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^ ! Abstract 

^: 
CI,: 

The properties of cosmic rays with energies above 10^ GeV have to be deduced from the spacetime 

. structure and particle content of the air showers which they initiate. In this review we summarize 

^ \ the phenomenology of these giant air showers. We describe the hadronic interaction models used to 

\ extrapolate results from collider data to ultra high energies, and discuss the prospects for insights 

\^ • into forward physics at the LHC. We also describe the main electromagnetic processes that govern 



the longitudinal shower evolution, as well as the lateral spread of particles. Armed with these two 
principal shower ingredients and motivation from the underlying physics, we provide an overview of 
some of the different methods proposed to distinguish primary species. The properties of neutrino 
interactions and the potential of forthcoming experiments to isolate deeply penetrating showers 
from baryonic cascades are also discussed. We finally venture into a terra incognita endowed with 
TeV-scale gravity and explore anomalous neutrino-induced showers. 
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I. INTRODUCTION 



A. Why another cosmic ray review? 

Some 40 years after the discovery of ultra high energy cosmic rays [1] , fundamental ques- 
tions regarding their origin and nature lack definitive answers. The highest primary energy 
measured thus far is ~ 10^^'^ GeV [2], corresponding to a nucleon-nucleon center-of-mass 
energy ^/s ~ 10^'^ GeV/\/~A, where A is the mass number of the primary particle. The 
existence of these particles, the most energetic known in the universe, challenges our current 
understanding of physics. From the perspective of astrophysics, one must identify where 
and how in the universe these particles obtain such high energies. A failure to uncover 
such mechanisms may lead one to postulate new physics to explain the phenomenon. From 
the perspective of particle physics, ultra high energy cosmic ray interactions are orders of 
magnitude beyond what can be achieved in current (and future) terrestrial collider exper- 
iments and may open a window to energy and kinematic regions previously unexplored in 
the study of fundamental interactions. From both perspectives, the tantalizing possibility 
of new physics that may be found in the study of ultra high energy cosmic rays continues 
to motivate current and future cosmic ray experiments. 

The literature abounds in reviews of experimental techniques for detection of cosmic ray 
air showers [3-9], as well as overviews of the physics of cosmic ray propagation and possible 
astrophysical and exotic origins [10-18]. This review follows a somewhat different path, 
focusing exclusively on cosmic ray phenomenology from the top of the atmosphere to the 
Earth's surface. The topics covered are viewed from the perspective of particle physics, and 
the reader is assumed only to possess a basic background in this field. We hope this article 
can provide a sort of bridge for high energy physicists interested in exploring some of the 
challenges facing upcoming ground- and space-based cosmic ray observatories. 



B. Cosmic ray observations 

For primary energy ^ 1 GeV the observed cosmic ray flux can be described by a series 
of power laws with the fiux falling about three orders of magnitude for each decade increase 
in energy. In the decade centered at \/s\^ ~ lO'^'"' QeV/\/~A, the spectrum steepens 
from to E~^'^, forming the feature commonly known as "the knee". The spectrum 

steepens further to E"^-^ above ~ 10^-^ GeV /\/A, and then flattens to E~'^-'^ at 

^/^l ankle ^ ^ GeV/vC4, formlng a feature known as "the ankle" [19]. Within the statistical 
uncertainty of existing data, which is large for E > 10^^ GeV, the tail of the spectrum is 
consistent with a simple extrapolation at that slope to the highest energies. Thus far, for 
Earth-based accelerators, the record holder for collisions with the highest energy per nucleon 
is the Tevatron, which countercirculates protons and antiprotons with ^/s ~ 1.8 TeV. This 
center-of-mass energy corresponds closely to that at the knee. The Large Hadron Collider 
(LHC), now under construction at CERN, will collide protons with protons at ^/s ~ 14 TeV. 
This impressive energy is still about a factor of 50 smaller than the center-of-mass energy 
of the highest energy cosmic ray so far observed, assuming A = 1. 

For primary cosmic ray energies above 10^ GeV, the flux becomes so low that direct 
detection of the primary using devices in or above the upper atmosphere is, for all practical 
purposes, impossible. Fortunately, in such cases the primary particle has enough energy to 
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initiate a particle cascade in the atmosphere large enough that the products are detectable 
at ground. There are several techniques which can be employed in detecting these extensive 
air showers (EAS), ranging from sampling of particles in the cascade to measurements of 
fluorescence, Cerenkov or radio emissions produced by the shower. 

The most commonly used detection method involves sampling the shower front at a 
given altitude using an array of sensors spread over a large area. Sensors, such as plastic 
scintillators or Cerenkov detectors are used to infer the particle density and the relative 
arrival times of the shower front at different locations; from this, one can estimate the energy 
and direction of the primary cosmic ray. The spacing between stations determines the energy 
threshold for a vertical shower. The muon content is usually sought either by exploiting the 
signal timing in the surface sensors or by employing dedicated detectors which are shielded 
from the electromagnetic shower component. Inferring the primary energy from energy 
deposits at the ground is not completely straightforward, and involves proper modeling of 
both the detector response and the physics of the first few cascade generations. This second 
point is particularly subtle and will be the main subject of Sec. 11. 

Another highly successful air shower detection method involves measurement of the lon- 
gitudinal development of the cascade by sensing the fluorescence hght produced via inter- 
actions of the charged particles in the atmosphere. As an extensive air shower develops, it 
dissipates much of its energy by exciting and ionizing air molecules along its path. Excited 
nitrogen molecules fluoresce producing radiation in the 300 - 400 nm ultraviolet range, to 
which the atmosphere is quite transparent. Under favorable atmospheric conditions EAS can 
be detected at distances as large as 20 km, about 2 attenuation lengths in a standard desert 
atmosphere at ground level, though observations can only be made on clear moonless nights, 
yielding a duty cycle of about 10%. The shower development appears as a rapidly moving 
spot of light whose angular motion depends on both the distance and the orientation of the 
shower axis. The integrated light signal is proportional to the total energy deposited in the 
atmosphere. Systematic errors can arise from a variety of sources, including uncertainties 
in the nitrogen fluorescence induced by the particle beam [20, 21], as well as uncertainties 
in the atmospheric conditions at the time the fluorescence measurements are taken [22]. 

The first measurements of ultra high energy cosmic rays were carried out by Linsley at Vol- 
cano Ranch (35°09'N, 106°47'W) in the late 1950's [23] using an array of scintillation coun- 
ters. More recent experiments using surface detection techniques include Haverah Park in 
England (53°58'N, 1°38'W) [24], Yakutsk in Russia (62°N, 130°E) [25], the Sydney University 
Giant Airshower Recorder (SUGAR) in Austraha (30°32'S, 149°43'E) [26], and the Akeno 
Giant Air Shower Array (AGASA), about 100 km west of Tokyo (38°4rN, 138°30'E) [27, 28]. 
The fluorescence method has been used by the Fly's Eye experiment [29, 30], as well as its 
up-scoped descendant High Resolution Fly's Eye experiment (HiRes) [31], operating at the 
Dugway proving ground in the Utah desert (40°N, 112°W). 

Over the next few years, the best observations of the extreme end of the cosmic ray 
spectrum will be made by the Pierre Auger Observatory (PAO) [32], which is currently 
operational in Malargiic, Argentina (35°12'S, 69°12'W) and is in the process of growing to 
its flnal size of 3000 km^. A twin site is pending for the Northern hemisphere, and together 
the two observatories will have an acceptance of 14000 km^ sr above 10^*^ GeV for zenith 
angles below 60°. The PAO works in a hybrid mode, and when complete will comprise 24 
fluorescence detectors overlooking a ground array of 1600 water Cerenkov detectors. During 
clear, dark nights, events are simultaneously observed by fluorescence light and particle 
detectors, allowing powerful reconstruction and cross-calibration techniques. Simultaneous 
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FIG. 1: Upper end of the cosmic ray energy spectrum as observed by AGASA [37] and 
HiRes [40] /Fly's Eye [41]. 

observations of showers using two distinct detector methods will also help to control the 
systematic errors that have plagued cosmic ray experiments to date. Moreover, each site 
of PAO reaches ~ 15 km^we sr of target mass around 10^° GeV [33], which is comparable 
to dedicated neutrino detectors being planned.^ This renders PAO a neutrino telescope 
operating in an energy regime complimentary to existing and upcoming facilities. The 
characteristics of neutrino-induced showers are discussed in Sec. V. 

Space-based experiments are also in the offing, the most thoroughly planned being the 
Extreme Universe Space Observatory (EUSO) [34, 35]. EUSO comprises a single fluorescence 
eye, and is scheduled to fly aboard the International Space Station for more than three years. 
After taking account of the 10% duty cycle, this experiment will image a vast volume of 
750 km^we sr. 

In recent years, a somewhat confused picture vis-a-vis the energy spectrum and ar- 
rival direction distribution has been emerging. Since 1998, the AGASA Collaboration has 
consistently reported [36, 37] a continuation of the spectrum beyond the expected Greisen- 



we = water equivalent. 
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Zatsepin-Kuzmin (GZK) cutoff [38, 39], wliicli sliould arise at about 10^°'^ GeV if cosmic 
ray sources are at cosmological distances. This theoretical feature of the spectrum is mainly 
a consequence of interactions of the primary cosmic ray with the microwave background 
radiation. In contrast, the most recent results from HiRes [40] describe a spectrum which is 
consistent with the expected GZK feature. The discrepancy between the 2 estimated fluxes 
is shown in Fig. 1. Several analyses were done in trying to understand this difference. In par- 
ticular, since the quoted systematic errors in the energy estimate are in the neighbourhood 
of 30%, it was argued [42] that if the AG AS A Collaboration overestimates their energies by 
30% and the HiRes Collaboration underestimates them by about the same amount, then 
the two spectra can be brought into reasonable agreement within statistical errors. 

Although there seems to be a remarkable agreement among experiments on predictions 
about isotropy on large scale structure [43, 44], this is certainly not the case when considering 
the two-point correlation function on a small angular scale. The AGASA Collaboration 
reports observations of event clusters which have a chance probability smaller than 1% to 
arise from a random distribution [45, 46]. Far from confirming the AGASA results, the 
recent analysis reported by the HiRes Collaboration showed that their data are consistent 
with no clustering among the highest energy events [47, 48]. The discovery of such clusters 
would be a tremendous breakthrough for the field, but the case for them is not yet proved. 
Special care must be taken when computing the statistical significance in such an analysis. 
In particular, it is important to define the search procedure a priori in order to ensure one 
does not inadvertently perform "trials" by studying the data before deciding upon the cuts. 
Very recently, with the aim of avoiding accidental bias on the number of trials performed 
in selecting the angular bin, the original claim of the AGASA Collaboration [45] was re- 
examined considering only those events observed after the original claim [49]. This study 
showed that the evidence for clustering in the AGASA data set is weaker than was previously 
supposed, and is consistent with the hypothesis of isotropically distributed arrival directions. 

The confusing experimental situation regarding the GZK feature and the small-scale 
clustering in the distribution of arrival direction should be resolved in the near future by the 
PAO, which will provide not only a data set of unprecedented size, but also the machinery 
for controlling some of the more problematic systematic uncertainties. As we will discuss in 
this review, however, the task of identifying primary species is more challenging. The rest 
of this article is organized as follows. In Sec. II, we describe the phenomenology of hadronic 
interactions with the goal of providing an overview of the main systematic uncertainties 
hindering the determination of the primary energy from observations at the ground. In the 
following section, we discuss the electromagnetic processes responsible for generating the 
great majority of particles in the shower. Armed with these two principal shower ingredients, 
we then discuss observables that are accessible to experiment. Next, in Sec. IV, we describe 
how these observables are used to infer the primary composition. In Sec. V we summarize 
properties of neutrino interactions and discuss observables in the deeply penetrating showers 
they produce. Since the expected rate of such events is very low, any enhancement of 
the cross section due to physics beyond the electroweak scale should be evident in this 
channel. At the end we allow ourselves to venture into speculative territory and discuss 
the experimental signatures of neutrino interactions in scenarios with TeV-scale gravity. 
Before getting underway, we first briefly summarize the main features of the atmospheric 
"calorimeter." 
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FIG. 2: Slant depths corresponding to various zenith angles 9 considering the curvature of the 
Earth. 

C. Nature's calorimeter 

Unlike man-made calorimeters, the atmosphere is a calorimeter whose properties vary in 
a predictable way with altitude, and in a relatively unpredictable way with time. Begin- 
ning with the easier of the two variations, we note that the density and pressure depend 
strongly on the height, while the temperature does not change by more than about 30% 
over the range 0-100 km above sea level. Therefore we can get a reasonable impression 
of the density variation by assuming an isothermal atmosphere, in which case the density 
Paxm{h) ^ PqC^^/^'^, where po ~ 1-225 kg/m^ and /iq = RT/(ng) fa 8.4 km is known as 
the scale-height of the atmosphere, R being the ideal gas constant, /i the average molecular 
weight of air, g the acceleration due to gravity and T ^ 288 K. Of course, reading out 
such a natural calorimeter is complicated by the effects of varying aerosol and molecular 
attenuation and scattering. 

The quantity that most intuitively describes the varying density of the atmospheric 
medium is the vertical atmospheric depth, Xy{h) — /^Patm(^) d^, where z is the height. 
However, the quantity most relevant in air shower simulations is the slant depth, X, which 
defines the actual amount of air traversed by the shower. The variation of the slant depth 
with zenith angle is shown in Fig. 2. If the Earth curvature is not taken into account, then 
X — Xy{h)/ cos6, where 6 is the zenith angle of the shower axis. For 9 < 80°, the error 
associated with this approximation is less than 4%. 

The atmospheric medium is endowed with a magnetic field. In general, the geomagnetic 
field is described by 3 parameters, its strength \B\, its inclination l, and its declination 5. 

— * 

The inchnation is defined as the angle between the local horizontal plane and the S-field. 
The declination is defined as the angle between the horizontal component of the field B± 
(i.e., perpendicular to the arrival direction of the air shower) and the geographical North 
(direction of the local meridian). The angle l is positive when B points downward and S is 
positive when B± is inclined towards the East. 
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II. HADRONIC PROCESSES 



Uncertainties in hadronic interactions at ultra high energies constitute one of the most 
problematic sources of systematic error in analysis of air showers. This section will explain 
the two principal schools of thought for extrapolating collider data to ultra high energies. We 
start with a general description of pp collisions within the eikonal model. Next, we consider 
the specific case of hadronic interactions in the atmosphere, and discuss the most widely 
used Monte Carlo codes. Finally, we study the potential of present and future accelerators to 
produce data valuable for understanding extensive air shower physics. Particular emphasis 
is placed on measurements of interaction processes at extreme forward directions and the 
cross sections for diffr active interactions. 

A. Low-Pj, jet physics beyond collider energies 

Soft multiparticle production with small transverse momenta with respect to the collision 
axis is a dominant feature of most hadronic events at center-of-mass energies 10 GeV < 

^/s < 50 GcV (see e.g., [50, 51]). Despite the fact that strict calculations based on ordinary 
QCD perturbation theory are not feasible, there are some phenomeno logical models that 
successfully take into account the main properties of the soft diffractive processes. These 
models, inspired by 1/A^ QCD expansion are also supplemented with generally accepted 
theoretical principles like duality, unitarity, Regge behavior, and parton structure. The 
interactions are no longer described by single particle exchange, but by highly complicated 
modes known as Reggeons. Up to about 50 GeV, the slow growth of the cross section with 
y/s is driven by a dominant contribution of a special Reggeon, the Pomeron. 

At higher energies, semihard interactions arising from the hard scattering of partons 
that carry only a very small fraction of the momenta of their parent hadrons can compete 
successfully with soft processes [52-59]. These semihard interactions lead to the minijct 
phenomenon, i.e., jets with transverse energy [Et = \Pj.\) much smaller than the total 
center-of-mass energy. Unlike soft processes, this low-p^ jet physics can be computed in 
perturbative QCD. The parton-parton minijet cross section is given by 

^QCD(.,pr') ^J2f^f^ d\i\ §f x,n{x,, x,f,{x,, , (1) 

■ J xi J X2 Jq2 . a\t\ 

tjj ^min I I 

where Xi and X2 are the fractions of the momenta of the parent hadrons carried by the 
partons which collide, daij/d\t\ is the cross section for scattering of partons of types i and j 
according to elementary QCD diagrams, fi and fj are parton distribution functions (pdf's), 
s — X1X2S and —i — s{l — cos'j9*)/2 = are the Mandelstam variables for this parton- 
parton process, and the sum is over all parton species. Here, 

p,^E^^ sin^jet = ^ sinr, (2) 

and 

P„=£;J,t COS^jet (3) 

where Ej^^ is the energy of the jet in the lab frame, i^jet the angle of the jet with respect to 
the beam direction in the lab frame, and -&* is the angle of the jet with respect to the beam 
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FIG. 3: Kinematic x-Q^ plane accessible to the HI and ZEUS experiments at HERA and the 
region accessible to fixed-target experiments. The inelasticity y = (1 — cost?*)/2 is also shown. 
This figure is courtesy of Max Klein. 

direction in the center-of-mass frame of the elastic parton-parton collision. This implies 
that for small '(9*, ~ . The integration limits satisfy Q^in < I'^l < ^I'^i with Qmin the 
minimal momentum transfer. 

A first source of uncertainty in modeling cosmic ray interactions at ultra high energy 
is encoded in the extrapolation of the measured parton densities several orders of magni- 
tude down to low X. Primary protons that impact on the upper atmosphere with energy 
> 10^"^ GeV, yield partons with x = 2p*/y/s > 1X1^,1^ ~ 10"'^, whereas current data on 

quark and gluon densities are only available for x > 10~^ to within an experimental accuracy 
of 3% for ^ 20 GeV^ [60]. In Fig. 3 we show the region of the x — plane probed by 
HI, ZEUS,^ and fixed target experiments. Moreover, application of HERA data to baryonic 
cosmic rays assumes universality of the pdf's. The QCD factorization conjecture, which 
is essentially equivalent to the Ingleman-Schlein model [61], posits that the parton-parton 
minijet cross section of Eq. (1) can always be written in a form which factorizes the par- 
ton densities and the hard interaction processes irrespective of the order in perturbation 



^ In the pe^ storage ring HERA at DESY, 27.6 GeV e^'s are collided on 820 GeV p's, and data are recorded 
by two experiments, HI and ZEUS. These collisions correspond to y/s w 300 GeV, or equivalently a lepton 
energy w 47 TeV in the proton rest frame. 
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theory and the particular hard process. This conjecture holds in the limit ^ Aqcd, 
where Aqcd ~ 200 MeV is the QCD renormalization scale. However, a severe breakdown 
of the factorization conjecture has been observed when using the pdf's obtained by the 
HERA experiments to predict diffractive jet production in hadron-hadron interactions at 
the Tevatron [62]. 

For large and not too small x, the Dokshitzer-Gribov-Lipatov-Altarelli-Parisi 
(DGLAP) equations [63-66] 



d\nQ^\g{x,Q^)J 27r P99J \9i^^Q 

successfully predict the dependence of the quark and gluon densities {q and g, respec- 
tively). Here, = Qg/lAir), with Qs the strong coupling constant. The splitting functions 
Pij indicate the probability of finding a daughter parton i in the parent parton j with a 
given fraction of parton j momentum. This probability will depend on the number of split- 
tings allowed in the approximation. In the double-leading-logarithmic approximation, limit 
[ln(l/a;), \n{Q'^ / Aq^^^)] 00, the DGLAP equations predict a steeply rising gluon density, 
xg ~ a;~°'^, which dominates the quark density at low x, in agreement with experimental 
results obtained with the HERA collider [67]. Specifically, as can be seen in Fig. 4, HERA 
data are found to be consistent with a power law, xg{x,Q^) ~ x~^^, with an exponent 
Ah between 0.3 and 0.4 [68]. However, it is easily seen using geometrical arguments that 
the rapid growth of the gluon density at low x would eventually require corrections to the 
evolution equations [69]. 

The high energy minijet cross section is then determined by the small-x behavior of the 
parton distributions or, rather, by that of the dominant gluon distribution (via the lower 
limits of the Xi, X2 integrations) 

ctqcb{s,pT') ^ d\i\ x,g{x,Ai\) x,g{x,Ai\) . (5) 

J Xi J X2 Jq2, d\t\ 

mm ' ' 

A naive estimate of the ctqcd behavior at high energies can be obtained via extrapolation 
of xg oc x~^^ to small x in Eq. (5). Within this approximation it is sufficiently accurate 
to keep only the leading contribution of the differential cross section for gg scattering {i.e., 
da/d\i\ oc |f|~^), and so Eq. (5) becomes [70] 

^qcd(s) oc ^ /' ^ ^ Hs/s,) , (6) 

J2Q2. u Xi J2Q2. u X2 

^mm' ^mm' 

where sq is a normalization constant. One caveat is that the inclusive QCD cross section 
given in Eq. (6) is a Born approximation, and therefore automatically violates unitarity. 

The procedure of calculating the inelastic cross section from inclusive cross sections is 
known as unitarization. In the eikonal model [71-74] of high energy hadron-hadron scat- 
tering, the unitarized (elastic, inelastic, and total) cross section, assuming a real eikonal 
function, is given by 



Cinel 



J d% [1 - exp [-x.„a(«, ^) - Xha.d(«, ^)] }' , (7) 
j d%[l- exp [-2x.„,, (s, h) - 2x,.,, (s, h)] } , (8) 
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FIG. 4: Gluon density for x > 10"^ and = 20 GeV^, as measured by HI, ZEUS and NMC 
experiments. 
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(9) 



where the scattering is compounded as a sum of QCD ladders via hard and soft processes 
through the eikonals Xhard X^oit ■ should be noted that we have ignored spin- dependent 
effects and the small real part of the scattering amplitude, both good approximations at 
high energies (see e.g., [75]). Now, if the eikonal function, x(s, b) = X^oai^' ^) + X^^rdi^' ^) ~ 
A/2, indicates the mean number of partonic interaction pairs at impact parameter b, the 
probability pn for having n independent partonic interactions using Poisson statistics reads. 



p„ = (AVn!) 



-A 3 



Therefore, the factor 1 — e = Yl'^=i Pn in Eq. (8) can be interpreted 



semiclassically as the probability that at least 1 of the 2 protons is broken up in a collision at 
impact parameter b. With this in mind, the inelastic cross section is simply the integral over 
all collision impact parameters of the probability of having at least 1 interaction, yielding a 
mean minijet multiplicity of (njet) ~ o'Qco/o'mei [78]. The leading contenders to approximate 
the (unknown) cross sections at cosmic ray energies, sibyll [79] and QGSJET [80], share the 



This relation can be derived within a field theoretical context [76] using the Abramovski-Gribov-Kancheli 
(AGK) cutting rules [77]. 
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eikonal approximation but differ in tfieir ansdtse for tlie eikonals. In botfi cases, tfie core of 
dominant scattering at very high energies is the parton-parton minijet cross section given 
in Eq. (1), 

X,..a-l<^QCB{s,pT^)A{s,b), (10) 

where the normahzed profile function, J dPh A{s, h) — 1, indicates the distribution of partons 
in the plane transverse to the collision axis. 

In the QGSJET-like models, the core of the hard eikonal is dressed with a soft-pomeron 
pre-evolution factor. This amounts to taking a parton distribution which is Gaussian in the 

— * 

transverse coordinate distance 

e-\i\VR\s) 

where -R^(s) ~ 4/?q + 4a^g ln^(s/so), with a'^f^ ^ 0.11 GcV"^. Fits to collider data have 
been carried out [81] using a Gaussian profile function with energy-independent width, Rq, 
which was allowed to vary in the fit. Under the assumption that the partons contributing 
to jet production are uniformly distributed in the transverse space all over the proton, 
one can obtain a reasonable fit to the data with Rq = 3.5 GeV^^ and p'^^"^ = 3.5 GeV. 
However, if one allows for the possibility of parton clustering, Rq shrinks to 1.5 GeV~^ with 
^cutoff _ 2 5 GeV. This leads to a smaller rise of the cross section with energy. In fact, 
the CDF Collaboration has reported [82] measurements which may indicate that partons 
are distributed in clusters inside the proton. Specifically, measurements of the 4-jet to 2-jet 
ratio for a jet transverse energy cutoff of 5 GeV when conveniently express in term of the 
effective cross section [83] lead to 



1 W2- 



■jetj 



12 



14.5± 1.7l^-^ mb. (12) 

^ <^4-jet 

Within the eikonal unitarization, this corresponds to cTcfr = SnR^. From Eq. (12), Rq !=s 
1.5 GeV~^, which is consistent with the clustering hypothesis. 

In SiBYLL-like models, the transverse density distribution is taken as the Fourier transform 
of the proton electric form factor, resulting in an energy-independent exponential (rather 
than Gaussian) fall-off of the parton density profile with The main characteristics of the 
pp cascade spectrum resulting from these choices are readily predictable: the harder form 
of the SIBYLL form factor allows a greater retention of energy by the leading particle, and 
hence less available for the ensuing shower. Consequently, on average SiBYLL-like models 
predict a smaller multiplicity than QGSJET-like models (see e.g. [84-87]). 

At high energy, Xsoh ^ Xhard' inelastic cross section is dominated by the 

hard eikonal. For impact parameters larger than some threshold, bg, where Xhardl-^' ^s) 3> 1, 
the damping from the exponential term in the Gaussian profile function of Eq. (11) is so 
strong that any increase in itqcd does not significantly alter the contribution to the inelastic 
cross section from the region where \b\ < bg. At high energy, with the appropriate choice 
of normalization, the cross section in Eq. (6) can be well-approximated by a power law. 
Hence, by taking uqcd ^ s^^, one fixes b^ ~ 4Q;gg Ah ln^(s/so) [86]. This implies that the 
growth of the inelastic cross section according to QGSJET-hke models is given by 

(Jinei ~ J d'^b e{bs - \b\) = nbl ~ Ana'^fi Ah ln^(s/so) ~ 0.52 Ah ln^(s/so) mb . (13) 
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FIG. 5: Energy dependence of the pp inelastic cross section as predicted by Eqs. (13) and (14) 
with 0.3 < Ah < 0.4. The darkly shaded region between the solid lines corresponds to the model 
with Gaussian parton distribution in b. The region between the dashed-dotted lines corresponds 
to the model with exponential fall-off of the parton density in b. In both cases the cross sections 
are normalized to reproduce the data (*) [89] from the CERN Intersecting Storage Ring (ISR) at 
30 GeV. Also shown are estimates [90] of the inelastic pp cross section as derived from measurements 
of the inelastic p-air cross section by the AGASA (■) [91] and the Fly's Eye (•) [92] experiments. 

For SIBYLL-Iike models, where Eq. (11) is replaced by the Fourier transform of the proton 
electric form factor, the growth of the inelastic cross section also saturates the In^ s Froissart 
bound [88], but with a multiplicative constant which is larger than the one in QGSJET-like 
models [86]. Namely, 



Figure 5 illustrates the large range of predictions for pp inelastic cross section which 
remain consistent with HERA data. When the two leading order approximations discussed 
above are extrapolated to higher energies, both are consistent with existing cosmic ray data. 
Note, however, that in both cases the range of allowed cross-sections at high energy varies 
by a factor of about 2 to 3. A point worth noting at this juncture: a number of approaches 
have been used to extract the pp cross section from cosmic ray shower data [93-96]. The 
points in Fig. 5 correspond to the most up-to-date estimate [90]. 

There are three event generators, sibyll [79], QGSJET [80], and dpmjet [97] which are 
tailored specifically for simulation of hadronic interactions up to the highest cosmic ray 




(14) 
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energies.^ The latest versions of these packages are SIBYLL 2.1 [99], QGSJET 01 [100], and 
DPMJET III [101]; respectively. In QGSJET, both the soft and hard processes are formulated 
in terms of Pomeron exchanges. To describe the minijets, the soft Pomeron mutates into a 
"semihard Pomeron" , an ordinary soft Pomeron with the middle piece replaced by a QCD 
parton ladder, as sketched in the previous paragraph. This is generally referred to as the 
"quasi-eikonal" model. In contrast, SIBYLL and dpmjet follow a "two channel" eikonal 
model, where the soft and the semi-hard regimes are demarcated by a sharp cut in the 
transverse momentum: SIBYLL 2.1 uses a cutoff parametrization inspired in the double 
leading logarithmic approximation of the DGLAP equations, 

= P° + 0.065 GeV exp[0.9 Vh[2] , (15) 
whereas dpmjet III uses an ad hoc parametrization for the transverse momentum cutoff 

(^) =pI + 0.12 GeV [logio(Vi/50GeV)]3 , (16) 

where p° = 2.5 GeV [68]. 

The transition process from asymptotically free partons to colour-neutral hadrons is de- 
scribed in all codes by string fragmentation models [102]. Different choices of fragmentation 
functions can lead to some differences in the hadron multiplicities. However, the main 
difference in the predictions of QGSJET-like and SIBYLL-Iike models arises from different 
assumptions in extrapolation of the parton distribution function to low energy. 



B. Hadronic interactions in the Earth's atmosphere 

Now we turn to nucleus-nucleus interactions, which cause additional headaches for event 
generators which must somehow extrapolate pp interactions in order to simulate the proton- 
air collisions of interest. All the event generators described above adopt the Glauber for- 
malism [71], which is equivalent to the eikonal approximation in nucleon-micleon scattering, 
except that the nucleon density functions of the target nucleus are folded with that of the 
nucleon. The inelastic and production cross sections read: 



Cinel 



"^prod 



exp 



<7tot TA{h) 



C^inel TA{h) 



(17) 



(18) 



where Ta(^) is the transverse density of hadronic matter of the target nucleus folded with 
that of the projectile hadron. Here, (Jinei and atot are given by Eqs. (8) and (9), respectively. 
The p-air inelastic cross section is the sum of the "quasi-elastic" cross section, which corre- 
sponds to cases where the target nucleus breaks up without production of any new particles, 
and the production cross section, in which at least one new particle is generated. Clearly 
the development of EAS is mainly sensitive to the production cross section. Overall, the 
geometrically large size of nitrogen and oxygen nuclei dominates the inclusive proton-target 



Additionally, a new event generator, nexus [98], is available for simulation of the region below ^/s 
10^ GeV. 
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FIG. 6: Left panel: The slowly rising curves indicate the mean inelasticity in proton air collisions 
as predicted by QGSJET and sibyll. The falling curves indicate the proton mean free path in the 
atmosphere. The data point is from Fly's Eye measurements [92]. Right panel: Mean multiplicity 
of charged secondary particles produced in inelastic proton-air collisions processed with QGSJET 
and SIBYLL. 

cross section, and as a result the disagreement from model-dependent extrapolation is not 
more than about 15%. 

The event generators also make different choices in their handling of nucleus-air col- 
lisions [103, 104]. Models of nucleus- nucleus interactions are particularly important to 
describe the first few generations of secondaries in cosmic ray showers produced by nu- 
clei. Measurements of proton-nucleus reactions at lower energies [105] suggested that the 
charged multiplicity from a soft production mechanism should simply scale with the number 
of nucleons that participate in the collision [106], thus allowing for comparison of different 
nuclear systems based on simple nucleon-nucleon superposition models. The particle den- 
sities are sensitive to the relative contributions of soft and hard processes [107, 108]. More 
recent experimental input suggests a simple superposition model is not completely realistic. 
Specifically, RHIC ^ data have shown that the observed central particle densities [109-111] 
are smaller than predictions from conventional multi-string models, with differences of 20% 
- 30% [112, 113]. To reduce the multiplicity in the models, the percolation process, which 
leads with increasing density to more and more fusion strings, has been proposed [114-116]. 
Very recently, the data from d-Au collisions collected by the PHOBOS Collaboration [117] 
were used to improve the event generator dpmjet III and bring the predicted multiplicity 
in line with the data [118]. 

So far the discussion has concerned p— air and nucleus-air interactions. Of course in air 
shower simulations, we are also concerned with vr-air interactions. Each approach discussed 



^ The Relativistic Heavy Ion Collider (RHIC) collides ultra relativistic ions at energies up to 0.2 TeV/N. 
RHIC has two large detectors, STAR and PHENIX and two smaller experiments: BRAHMS and PHO- 
BOS. 
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above handles up collisions using the same interaction model it uses for pp collisions. For 
energies of interest, both models predict, on average, a 7rp inelastic cross section about 20% 
smaller than the pp cross section [86] . 

Since the codes described above are still being refined, the disparity between them can 
vary even from version to version. At the end of the day, however, the relevant parameters 
boil down to two: the mean free path, A = (Jprod^)"^, and the inelasticity, K = 1 — 
-£'iead/-E'proj) whcrc Ti is the number density of atmospheric target nucleons, £^iead is the 
energy of the most energetic hadron with a long lifetime, and i^proj is the energy of the 
projectile particle. The first parameter characterizes the frequency of interactions, whereas 
the second one quantifies the energy lost per collision. Overall, SIBYLL has a shorter mean 
free path and a smaller inelasticity than QGSJET, as indicated in Fig. 6. Since a shorter 
mean free path tends to compensate a smaller inelasticity, the two codes generate similar 
predictions for an air shower which has lived through several generations. The different 
predictions for the mean charged particle multiplicity in proton-air collisions are shown in 
Fig. 6. Both models predict the same multiplicity below about 10^ GeV, but the predictions 
diverge above that energy. Such a divergence readily increases with rising energy. While 
QGSJET predicts a power law-like increase of the number of secondaries up to the highest 
energy, SIBYLL multiphcity exhibits a logarithmic growth. As it is extremely difficult to 
observe the first interactions experimentally, it is not straightforward to determine which 
model is closer to reality. In Sec. Ill, however, we will discuss observables which may offer 
a hint of which model better predicts overall shower characteristics. 

C. Measurements of forward processes at the LHC 

Interpretation of cosmic ray data suffers from the lack of knowledge of high energy 
hadronic interaction models. Hard interactions with high momentum transfer are calculable 
in perturbation theory using QCD. At present, collider experiments have mostly concen- 
trated on these hard processes in the central region, thereby excluding soft processes in the 
far-forward direction. These low momentum transfer processes, which are of great inter- 
est in the development of cosmic ray EAS, are not calculable from the fundamental QCD 
Lagrangian. 

Some guidance towards understanding hadronic processes in the forward direction may 
come directly from measurements of hadrons in airshowers [119]. However, the most useful 
experimental input in the forseeable future will likely come from the LHC. This machine, 
expected to become operational in 2007 or so, will provide pp collisions with ^/s = 14 TeV 
and luminosity up to L fti lO^'' cm~^ s~^ [120], as well as, a few years later, lead-lead 
ion collisions with ^/s — 1000 TeV. Two general-purpose experiments, ATLAS and CMS, 
presently cover up to 1 77] < 5. A dedicated heavy ion detector, ALICE, will also operate at 
this collider.^ 

The interesting low momentum transfer processes tend to populate the region at very 
small angles 1!} with respect to the beam direction. The distribution of pseudorapidity, 
77 = — lntan(t?/2), and the energy flow distribution are shown in Fig. 7. While the particle 
multiplicity is greatest in the low I77I region, it is clearly seen that the energy flow is peaked 



^ A 6-physics experiment, LHC-6, is also under construction at the LHC and will offer particle identification 
in the range 1.9 < 77 < 4.9. 
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at small production angles (large I77I). 

A study of diffraction, in pp as well as heavy ion collisions, must use detectors with ex- 
cellent forward acceptance to allow for a comparison with cosmic ray data. Dedicated runs 
of the LHC with lower luminosity (L = 10^^cm~^s~^) and specially tuned beam optics are 
planned to study these diffractive events. At present the only approved experiment at LHC 
with a capability of measuring, to some extent, very forward particles is TOTEM [121-123], 
which will comprise Roman pots placed on each side of the CMS interaction region and for- 
ward trackers which cover the pseudorapidity range 3.0 < rj < 6.8. It should be mentioned, 
however, that the fragmentation region that plays a crucial role in the development of EAS 
corresponds to pseudorapidity range 6 < |r;| < 10. 

The main goal of TOTEM is the measurement of elastic and total cross-sections with an 
expected precision of about 1%, in a luminosity independent manner. To calculate the total 
cross section in terms of the number of elastic and inelastic events measured by TOTEM, 
we can resort to the well-known optical theorem 

(Jtot = »m[/(0)] , (19) 

where f{i}) satisfies 

%i = i^^ = -|/(tf)P. (20) 

at s a\l s 

with i9 the angle of the scattered proton with respect to the beam direction. Squaring 
Eq.(19) we obtain 



.2 
tot 



167r W[/(0)] rfcTel 



3fJe2[/(0)] + ^Jm2[/(0)] dt 



IGtt [ dNei/dt]t=o ^^^^ 



t=0 



1 + p 



where C is the integrated luminosity. Now, following [124, 125], we can obtain the total 
cross section independently from C, by using utot = (cei + Cinei) = (A^ei + Mnei)/>C, 

_ IGtt [dNjdt]t=o 

Here, A^ei and A^inei are the numbers of elastic and inelastic events, and p ~ 0.10 ±0.01 is the 
ratio between the real and imaginary parts of the forward scattering amplitude [126].^ The 
difficult aspect of this measurement is obtaining a good extrapolation of the cross section 
for low momentum transfer. Recall that —t = s{l — cos't?)/2 ~ s'd'^/i, and so t — > implies 
a measurement in the extreme forward direction. The TOTEM experiment aims to measure 
down to values of \t\ xlO~^ GeV^, which corresponds to ■»? 4.5 //rad [123]. The design 
for the pseudorapidity range 5.5 < I77I < 6.8 is under discussion in a joint CMS/TOTEM 
working group. A tungsten Cerenkov calorimeter known as CASTOR has also been proposed 
which would compliment the measurements of TOTEM and CMS for |?7| < 6.8 and facilitate 
simultaneous measurements of particle fiow in diffractive and non-diffractive events. 

The ATLAS experiment is planning to implement additional detectors to cover the for- 
ward diffractive regions with tracking and/or calorimetry [127], with proposed coverage for 
the region |t| pa 6 x 10""^ GeV 2. 



Note that the quoted value of p is an extrapolation to ^/s = 14 TeV, and may be measured by the LHC 
experiments. Otherwise, it will contribute to the uncertainty in atot- 



17 



6 
5 
4 
3 
2 
1 



.4 
.2 





LHC, inelastic collisions 

TOTEM, CMS 

< -. . > 

Tracking and Calorimetry 



0.8 
0.6 
0.4 
0.2 




ATLAS, CMS 
Calorimetry 



ATLAS, CM S 
Tracking 

-10 -7.5 -5 -2.5 2.5 



7.5 10 



FIG. 7: Pscudorapidity distributions of charged particles (upper panel) and of the energy flow 
(lower panel) for pp collisions at LHC [121]. 



In summary, existing event generators rely on theoretical extrapolations of existing data 
up to the energies near the GZK energy. There is a general consensus in the community that 
in order to understand the development of EAS at these extreme energies, new inputs from 
accelerator experiments are needed. A series of workshops have been organized to discuss 
what experimental inputs are most needed [128-130]; a preliminary list of the requirements 
includes [131]: (i) measurements of total and inelastic cross sections for pp, pA, AA (ii) 
measurements of the ratio between soft diffractive and semi-hard processes, (Jdiff/cinei, (Hi) 
measurement of inclusive final state hadrons in the two momentum ranges, 0.8 < x < 1.0 
and 0.1 < x < 0.8. 



III. ELECTROMAGNETIC PROCESSES 

In this section we describe the electromagnetic interactions of relevance in ultra high 
energy shower development. The most important processes are electron and muon 
bremsstrahlung and pair production. We also discuss the Landau- Pomeranchuk-Migdal 
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(LPM) effect, wliicli suppresses tlie cross sections for pair production and brcmsstralilung 
above roughly 10^° GeV, and plioton conversion in tfie geomagnetic field, which to a large 
degree compensates for the LPM effect in terms of shower observables. We comment further 
on shower observables such as the age parameter, Moliere radius and shower size within the 
context of the Nishimura, Kamata and Greisen (NKG) model. Finally, we discuss the exten- 
sion of the NKG formalism to the corresponding shower parameters describing the lateral 
spread and the longitudinal development of EAS initiated by hadrons. 

A. The electromagnetic component 

The evolution of an extensive air shower is dominated by electromagnetic processes. The 
interaction of a baryonic cosmic ray with an air nucleus high in the atmosphere leads to 
a cascade of secondary mesons and nucleons. The first few generations of charged pions 
interact again, producing a hadronic core, which continues to feed the electromagnetic and 
muonic components of the showers. Up to about 50 km above sea level, the density of 
atmospheric target nucleons is n ~ 10^° cm~^, and so even for relatively low energies, say 
i?7r± ~ 1 TeV, the probability of decay before interaction falls below 10%. Ultimately, the 
electromagnetic cascade dissipates around 90% of the primary particle's energy, and hence 
the total number of electromagnetic particles is very nearly proportional to the shower 
energy [132]. 

By the time a vertically incident 10^^ GeV proton shower reaches the ground, there 
are about 10^^ secondaries with energy above 90 kcV in the the annular region extending 
8 m to 8 km from the shower core. Of these, 99% are photons, electrons, and positrons, 
with a typical ratio of 7 to e+e" of 9 to 1. Their mean energy is around 10 MeV and 
they transport 85% of the total energy at ground level. Of course, photon-induced showers 
are even more dominated by the electromagnetic channel, as the only significant muon 
generation mechanism in this case is the decay of charged pions and kaons produced in 7-air 
interactions [133]. 

It is worth mentioning that these figures dramatically change for the case of very inclined 
showers. For a primary zenith angle, 9 > 70°, the electromagnetic component becomes 
attenuated exponentially with atmospheric depth, being almost completely absorbed at 
ground level. We remind the reader that the vertical atmosphere is ~ 1000 g/cm^, and is 
about 36 times deeper for completely horizontal showers (see Fig. 2). As a result, most of 
the energy at ground level from an inclined shower is carried by muons. 

In contrast to hadronic collisions, the electromagnetic interactions of shower particles can 
be calculated very accurately from quantum electrodynamics. Electromagnetic interactions 
are thus not a major source of systematic errors in shower simulations. The first compre- 
hensive treatment of electromagnetic showers was elaborated by Rossi and Greissen [134]. 
This treatment was recently cast in a more pedagogical form by Gaisser [135], which we 
summarize in the subsequent paragraphs. 

The generation of the electromagnetic component is driven by electron bremsstrahlung 
and pair production [136]. Eventually the average energy per particle drops below a critical 
energy, eo, at which point ionization takes over from bremsstrahlung and pair production 
as the dominant energy loss mechanism. The energy loss rate due to bremsstrahlung 
radiation is nearly proportional to their energy, whereas the ionization loss rate varies only 
logarithmically with the energy. Though several different definitions of the critical energy 
appear in the literature [137] , throughout this review we take the critical energy to be that 
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at which the ionization loss per radiation lenght is equal to the electron energy, yielding 
eo = 710 MeV/(Zcff + 0.92) ~ 86 MeV [138].^ The changeover from radiation losses to 
ionization losses depopulates the shower. One can thus categorize the shower development 
in three phases: the growth phase, in which all the particles have energy > eo; the shower 
maximum, X^g^^; and the shower tail, where the particles only lose energy, get absorbed or 
decay. 

The relevant quantities participating in the development of the electromagnetic cascade 
are the probability for an electron of energy E to radiate a photon of energy k = yE and 
the probability for a photon to produce a pair e'^e" in which one of the particles (hereafter 
e~) has energy E = xk. These probabilities are determined by the properties of the air and 
the cross sections of the two processes. 

In the energy range of interest, the impact parameter of the electron or photon is 
larger than an atomic radius, so the nuclear field is screened by its electron cloud. In 
the case of complete screening, where the momentum transfer is small, the cross section for 
bremsstrahlung can be approximated by [140] 

' - ' o-o?/ + 2/ ' (23) 



dk XoNAk V3 3 

where A^s is the effective mass number of the air, Xq is a constant, and Na is Avogadro's 
number. In the infrared limit {i.e., y <^ 1) this approximation is inaccurate at the level of 
about 2.5%, which is small compared to typical experimental errors associated with cosmic 
air shower detectors. Of course, the approximation fails as y ^ 1, when nuclear screening 
becomes incomplete, and as y — > 0, at which point the LPM and dielectric suppression 
effects become important, as we discuss below. 

Using similar approximations, the cross section for pair production can be written as [140] 

— pa 1 X + -X . (24) 

dE XqNa V 3 3 y ^ ' 

The similarities between this expression and Eq. (23) are to be expected, as the Feynman 
diagrams for pair production and bremsstrahlung are variants of one another. 

The probability for an electron to radiate a photon with energy in the range {k, k + dk) 
in traversing dt — dX/Xo of atmosphere is 

dae^y XqNa , 4 1 - y 



dkdt^{y+ -]dydt, (25) 

dk A,s V ^ y J ^ ^ 

whereas the corresponding probability density for a photon producing a pair, with electron 
energy in the range {E, E + dE), is 



da^^e+e- XqNa 



dE 



ieflf 



dEdt^ (l--x+-x^) dx dt . (26) 



For altitudes up to 90 km above sea level, the air is a mixture of 78.09% of N2, 20.95% of O2, and 0.96% of 
other gases [139]. Such a mixture is generally modeled as an homogeneous substance with atomic charge 
and mass numbers Z^g = 7.3 and A^s = 14.6, respectively. 
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The total probability for pair production per unit of Xq follows from integration of Eq. (26), 

As can be seen from Eq. (25), the total probability for bremsstrahlung radiation is log- 
arithmically divergent. However, this infrared divergence is eliminated by the interference 
of bremsstrahlung amphtudes from multiple scattering centers. This collective effect of the 
electric potential of several atoms is known as the Landau-Pomeranchuk-Migdal (LPM) ef- 
fect [141, 142]. Of course, the LPM suppression of the cross section results in an effective 
increase of the mean free path of electrons and photons. This effectively retards the de- 
velopment of the electromagnetic component of the shower. It is natural to introduce an 
energy scale, -Elpm, at which the inelasticity is low enough that the LPM effect becomes 
significant [143]. Below S'lpm, the energy loss rate due to bremsstrahlung is roughly 

With this in mind, we now identify the constant Xq pa 36.7 g cm~^ with the radiation 
length in air, defined as the mean distance over which a high-energy electron loses 1/e of 
its energy, or equivalently 7/9 of the mean free path for pair production by a high-energy 

photon [137]. 

The experimental confirmation of the LPM effect at Stanford Linear Accelerator Center 
(SLAC) [144, 145] has motivated new analyses of its consequences in cosmic ray physics [146- 
150]. The most evident signatures of the LPM effect on shower development are a shift in the 
position of the shower maximum X^ax and larger fluctuations in the shower development. 

When considering the LPM effect in the development of air showers produced by ultra 
high energy cosmic rays, one has to keep in mind that the suppression in the cross sec- 
tions is strongly dependent on the atmospheric depth. ^ Since the upper atmosphere is very 
thin the LPM effect becomes noticeable only for photons and electrons with energies above 
EhPM ~ 10^° GeV. For baryonic primaries the LPM effect does not become important un- 
til the primary energy exceeds lO^^GeV. This is because the electromagnetic shower does 
not commence until after a significant fraction of the primary energy has been dissipated 
through hadronic interactions. To give a visual impression of how the LPM effect slows 
down the initial growth of high energy photon-induced showers, we show the average longi- 
tudinal shower development of 10^° GeV proton and 7-ray showers (generated using aires 
2.6.0 [151]) with and without the LPM effect in Fig. 8. 

At energies at which the LPM effect is important {viz., E > -Blpm), 7-ray showers will 
have already commenced in the geomagnetic field at almost all latitudes. This reduces the 
energies of the primaries that reach the atmosphere, and thereby compensates the tendency 
of the LPM effect to retard the shower development. The first description of photon inter- 
actions in the geomagnetic field dates back at least as far as 1966 [152], with a punctuated 
revival of activity in the early 1980's [153]. More recently, a rekindling of interest in the topic 
has led to refined calculations [154-156]. Primary photons with energies above 10^° GeV 
convert into e+e~ pairs, which in turn emit synchrotron photons. Regardless of the primary 



The same occurs for dielectric suppression, although the influence is not as important as for the LPM 
effect [148]. 
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FIG. 8: Average longitudinal shower developments of 10^^ GeV proton (dashed-dotted line) and 
7-rays with and without the LPM effect (solid and dotted lines, respectively). The primary zenith 
angle was set to 6 = 60°. The shadow area represents the intrinsic fluctuations of the showers. 
Larger fluctuations can be observed for 7-ray showers with the LPM effect, as expected. 

energy, the spectrum of the resulting photon "preshower" entering the upper atmosphere 
extends over several decades below the primary photon energy, and is peaked at energies 
below 10^° GeV [154]. The geomagnetic cooling thus switches on at about the same energy 
at which the LPM effect does, and thereby preempts the LPM-related observables which 
would otherwise be evident. Recent simulations [157] which include photon preshowering 
indicate that above ~ 10^^ GeV this effect tends to accelerate the shower development, 
shifting the X^ax to a smaller value than previous calculation suggested [158, 159], and into 
a range consistent with Xmax typical of a proton primary. 

The relevant parameter to determine both conversion probability and synchrotron emis- 
sion is E X B^, where E is the 7-ray energy and 5^ the transverse magnetic field. This 
leads to a large directional and geographical dependence of shower observables. Thus, each 
experiment has its own preferred direction for identifying primary gamma rays. For instance. 
Fig. 9 shows a map of the photon conversion probability in the geomagnetic field for all in- 
cident directions evaluated at the location of the HiRes experiment (|-B| = 0.53 G, l = 25°, 
and S = 14°) [155]. The smallest probabilities for conversion are found, not surprisingly, 
around the direction parallel to the local geomagnetic field. Note that this conversion-free 
region shrinks rapidly with increasing primary energy. A similar evaluation for the Southern 
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FIG. 9: Maps of gamma ray conversion probability in the geomagnetic field for several primary 
energies. Azimuths are as labeled, "N" denotes true north. The inner circles correspond to zenith 
angles 30°, 60° and horizon. Dashed curves indicate the opening angles of 30°, 60° and 90° to the 
local magnetic field [155]. 



Site of the Pierre Auger Observatory 
in [154]. 



\B\ 



0.25 G, L = -35°, and 5 = 86°) can be found 



1. Paper- and-pencil air shower modeling 

Most of the general features of an electromagnetic cascade can be understood in terms 
of the toy model due to Heitler [160]. In this model, the shower is imagined to develop 
exclusively via bremsstrahlung and pair production, each of which results in the conversion of 
one particle into two. As was previously discussed, these physical processes are characterized 
by an interaction length Xq. One can thus imagine the shower as a particle tree with branches 
that bifurcate every Xq, until they fall below a critical energy, eo, at which point energy 
loss processes dominate. Up to eo, the number of particles grows geometrically, so that after 
n — X/Xq branchings, the total number of particles in the shower is A?^ ~ 2". At the depth 
of shower maximum X^ax) all particles are at the critical energy, Cq, and the energy of the 
primary particle, Eq, is split among all the N-^^ = Eq/cq particles. Putting this together, 
we get: 

ln(£^o/eo) 



Y ~ 

^max ~ 



In 2 



(29) 



In real life, the combination of the LPM and geomagnetic effects introduces large fluctuations 
in the value of Xmax for photon showers. The prediction of this toy model roughly lies within 
the range of these fluctuations. 

Even baryon-induced showers are dominated by electromagnetic processes, so this toy 
model is still enlightening for such cases. In particular, for proton showers, Eq. (29) tells us 
that the Xmax scales logarithmically with primary energy, while N^^ax scales linearly. More- 
over, to extend this discussion to heavy nuclei, we can apply the superposition principle 
as a reasonable first approximation. In this approximation, we pretend that the nucleus 
comprises unbound nucleons, such that the point of first interaction of one nucleon is inde- 
pendent of all the others. Specifically, a shower produced by a nucleus with energy E^ and 
mass A is modeled by a collection of A proton showers, each with A~^ of the nucleus energy. 
Modifying Eq. (29) accordingly one easily obtains Xmax oc \n{Eo/A). 

While the Heitler model is very useful for imparting a first intuition regarding global 
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shower properties, the details of shower evolution are far too complex to be fully described 
by a simple analytical model. Full Monte Carlo simulation of interaction and transport 
of each individual particle is required for precise modeling of the shower development. At 
present two Monte Carlo packages are available to simulate EAS: CORSIKA (COsmic Ray 
Simulation for KAscade) [161] and aires (AIR shower Extended Simulation) [151]. Both 
programs provide fully 4-dimensional simulations of the air showers initiated by protons, 
photons, and nuclei. To simulate hadronic physics the programs make use of the event 
generators described in Sec. II. Propagation of particles takes into account the Earth's 
curvature and geomagnetic field. For further details on these codes the reader is referred 
to [162]. 

As a bridge between the first order approximation just described and a full-blown Monte 
Carlo treatment of air shower cascades, a hybrid method has recently been presented [163]. 
The approach is as follows. The first few interactions are treated using Monte Carlo event 
generators. The second step in the approach utilizes one-dimensional cascade equations up 
to the point where the lateral spread of the particles becomes non-negligible, then the output 
of the cascade equations is treated again with Monte Carlo. The method shows a reasonable 
agreement when compared with results of the two detailed simulation packages [164]. 



2. Electron lateral distribution function 

The transverse development of electromagnetic showers is dominated by Coulomb scatter- 
ing of charged particles off the nuclei in the atmosphere. The lateral development in electro- 
magnetic cascades in different materials scales well with the Moliere radius tm = EgXo/eo, 
which varies inversely with the density of the medium, 

X Patm(/iOL) 9.0 g/cm^ 

ru = ruihoh) tt^ ^ tt^ , (30) 

where Eg ^ 21 MeV [137] and the subscript OL indicates a quantity taken at a given 
observation level. 

Approximate calculations of cascade equations in three dimensions to derive the lateral 
structure function for a pure electromagnetic cascade in vertical showers were obtained by 
Nishimura and Kamata [165], and later worked out by Greisen [166] in the well-known NKG 
formula, 

where is the total number of electrons, r is the distance from the shower axis, and 

27rr(s,KG)r(4.5-2s,,J " ^ ^ 



For a primary of energy Eq, the so-called "age parameter", 

-1 

(33) 



^ ^ 2 ln(Eo/eo) 



t 



characterizes the stage of the shower development in terms of the depth of the shower in 
radiation lengths, i.e., t = Pa.tm{z) dz/X^. 
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FIG. 10: Electron (top) and muon (bottom) lateral distributions of a 10^*^ GeV vertical proton 
shower at different atmospheric altitudes. The solid lines are fits to the data using the NKG-like 
parametrization. The error bars are, in all cases, smaller than the points [167]. 

The NKG formula may also be extended to describe showers initiated by baryons [167]. 
In such an extension, one finds a deviation of behavior of the Moliere radius described 
in Eq. (30) when using a value of the age parameter which is derived from theoretical 
predictions for pure electromagnetic cascades. The need for a different age parameter to 
reproduce the electromagnetic component of hadronic induced showers has been addressed 
experimentally [168-175] and extensively studied by several authors [175-177]. It is possible 
to generalize the NKG formula for the electromagnetic component of baryon-induced show- 
ers by modifying the exponents in Eq. (31). From simulations, fits to lateral distribution 
functions (LDF) of electrons and positrons as a function of depth, t, yield an age parameter 
given by 

^ = 3(1 + ^)", (34) 

where the floating parameter j3 takes into account the above mentioned deviations from the 
theoretical value Sj^j^^ [167]. 

The modified NKG formula provides a good description of the e~^e~ lateral distribution 
at all stages of shower development for values of r sufficiently far from the hadronic core. 
Fortunately, this is the experimentally interesting region, since typical ground arrays can 
only measure densities at r > 100 m from the shower axis, where detectors are not saturated. 
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To illustrate the validity of this parametrization, we show in Fig. 10 (top) the Monte Carlo 
e^e~ density distributions corresponding to a single 10^° GeV proton shower at selected 
atmospheric depths. The total number of electrons obtained from the fit to each single 
shower is slightly lower than the true value due to the invalidity of the parametrization 
close to the shower core. It should be mentioned that an NKG-like formula can be used to 
parametrize the total particle's density observed in baryon-induced showers [178]. 

In the case of inclined showers, one normally analyzes particle densities in the plane 
perpendicular to the shower axis. Simply projecting distributions measured at the ground 
into this plane is a reasonable approach for near- vertical showers, but is not sufficient for 
inclined showers. In the latter case, additionally asymmetry is introduced because of both 
unequal attenuation of the electromagnetic components arriving at the ground earlier than 
and later than the core [167], and geometrical effects which also reduce the early compared 
to the late flux [179]. Moreover, deflections on the geomagnetic held become important for 
showers inclined by more than about 70°. 

In the framework of cascade theory, any effect coming from the influence of the atmosphere 
should be accounted as a function of the slant depth t [165]. Following this idea, a LDF 
valid at all zenith angles 9 < 70° can be determined by considering 

i'(^,C) = i sec^(l + /C cosC)"S (35) 

where C is the azimuthal angle in the shower plane, K, — KLq tan^, and /Co is a constant 
extracted from the fit [167, 180]. Then, the particle lateral distributions for inclined showers 
p(r, t') are given by the corresponding vertical LDF p(r, t) but evaluated at slant depth 
t'{9, Q where the dependence on the azimuthal angle is evident. 

For zenith angles d > 70°, the surviving electromagnetic component at ground is mainly 
due to muon decay and, to a much smaller extent, hadronic interactions, pair production and 
bremsstrahlung. As a result the lateral distribution follows that of the muon rather closely. 
In Fig. 11 the longitudinal development of the muon and electron components arc shown. 
It is evident from the figure that for very inclined showers the electromagnetic development 
is due mostly to muon decay [181, 182]. 

The consequences of the LPM effect and pair production in the geomagnetic field on the 
longitudinal cascade distribution initiated by photons were already discussed in this section. 
Since the lateral distribution of particles is strongly correlated with the development of the 
shower in the atmosphere, the LPM effect has consequences for the observed LDF at ground 
level. In particular, unconverted photons result in large fluctuations and steeper lateral 
profiles than nuclei [154]. 

In summary, the growth of the electromagnetic cascade is governed by bremsstrahlung 
and pair production. The mean free path for interactions via these processes depends on 
energy and atmospheric depth. Below 10^° GeV, each particle sees screened nuclei, while 
at higher energy collective effects suppress the cross section. On top of that, ultra high 
energy gamma ray interactions in the geomagnetic field also come into play, reducing the 
importance of the LPM cross section suppression. 

The lateral distribution of the electromagnetic component of a shower can be effectively 
parametrized. The well-known NKG lateral distribution function, which strictly applies 
to only purely electromagnetic showers, can be extended to describe not only the electro- 
magnetic portion of baryon-induced showers but also the signal produced by all particles 
reaching ground level. This provides a handle on one of the most useful shower observables 
available to surface arrays. 
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FIG. 11: Longitudinal development of muons and electrons as a function of the slant depth for 
10^^ GeV proton-induced showers. 

B. The muon component 

The muonic component of EAS differs from the electromagnetic component for two main 
reasons. First, muons are generated through the decay of cooled {Et,± < 1 TeV) charged 
pions, and thus the muon content is sensitive to the initial baryonic content of the primary 
particle. Furthermore, since there is no "muonic cascade", the number of muons reaching 
the ground is much smaller than the number of electrons. Specifically, there are about 
5 X 10*^ muons above 10 MeV at ground level for a vertical 10^^ GeV proton induced shower. 
Second, the muon has a much smaller cross section for radiation and pair production than 
the electron, and so the muonic component of EAS develops differently than does the elec- 
tromagnetic component. The smaller multiple scattering suffered by muons leads to earlier 
arrival times at the ground for muons than for the electromagnetic component. 

The ratio of electrons to muons depends strongly on the distance from the core; for 
example, the e~^e~ to fi'^fi~ ratio for a 10^^ GeV vertical proton shower varies from 17 to 1 at 
200 m from the core to 1 to 1 at 2000 m. The ratio between the electromagnetic and muonic 
shower components behaves somewhat differently in the case of inclined showers. For zenith 
angles greater than 60°, the e~ / fi'^ fi~ ratio remains roughly constant at a given distance 
from the core. As the zenith angle grows beyond 60°, this ratio decreases, until at ^ = 75°, 
it is 400 times smaller than for a vertical shower. Another difference between inclined 
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and vertical showers is that the average muon energy at ground changes dramatically. For 
horizontal showers, the lower energy muons are filtered out by a combination of energy loss 
mechanisms and the finite muon lifetime: for vertical showers, the average muon energy is 

1 GeV, while for horizontal showers it is about 2 orders of magnitude greater. The muon 
densities obtained in shower simulations using sibyll 2.1 fall more rapidly with lateral 
distance to the shower core than those obtained using QGSJET 01. This can be understood 
as a manifestation of the enhanced leading particle effect in sibyll, which can be traced to 
the relative hardness of the electromagnetic form factor profile function. The curvature of the 
distribution {d^Pn/ dr^) is measurably different in the two cases, and, with sufficient statistics, 
could possibly serve as a discriminator between hadronic interaction models, provided the 
primary species can be determined from some independent obscrvable(s) [183]. 

High energy muons lose energy through e+e~ pair production, muon-nucleus interaction, 
bremsstrahlung, and knock-on electron (5-ray) production [184]. The first three processes 
are discrete in the sense that they are characterized by high inelasticity and a large mean 
free path. On the other hand, because of its short mean free path and its small inelas- 
ticity, knock-on electron production can be considered a continuous process. The muon 
bremsstrahlung cross section is suppressed by a factor of {nie/m^y with respect to electron 
bremsstrahlung, see Eq. (23). Since the radiation length for air is about 36.7 g/cm^, and 
the vertical atmospheric depth is 1000 g/cm^, muon bremsstrahlung is of negligible impor- 
tance for vertical air shower development. Energy loss due to muon-nucleus interactions is 
somewhat smaller than muon bremsstrahlung. As can be seen in Fig. 12, energy loss by pair 
production is slightly more important than bremsstrahlung at about 1 GeV, and becomes 
increasingly dominant with energy. Finally, knock-on electrons have a very small mean free 
path (see Fig. 12), but also a very small inelasticity, so that this contribution to the energy 
loss is comparable to that from the hard processes. 

In addition to muon production through charged pion decay, photons can directly generate 
muon pairs, or produce hadron pairs which in turn decay to muons. In the case of direct 
pair production, the large muon mass leads to a higher threshold for this process than for 
electron pair production. Furthermore, QED predicts that fi^fi^ production is suppressed 
by a factor (me/mij)'^ compared the Bethe-Heitler cross section. The cross section for hadron 
production by photons is much less certain, since it involves the hadronic structure of the 
photon. This has been measured at HERA for photon energies corresponding to E'lab = 

2 X 10^ GeV [185, 186]. This energy is still well below the energies of the highest energy 
cosmic rays, but nonetheless, these data do constrain the extrapolation of the cross sections 
to high energies. To give an idea of the rates, at 100 GeV the cross section for 7 e^e~ is 
Pi 650 mb, i.e. much larger than the cross sections for hadronic interaction (f« 1.4 mb) or 
for muon pair production (r^ 0.015 mb). 

1. Muon lateral distribution function 

Now we consider the lateral distribution of the muon component of an extensive air 
shower. Unlike the electrons and photons, muons are relatively unaffected by multiple 
Coulomb scattering, and so their lateral distribution function (LDF) retains some informa- 
tion about the parent pion trajectories. In what follows, we first discuss a parameterization 
characterizing the muon LDF which is motivated only by the muon genealogy. After that, 
we add to the discussion the effect of the geomagnetic field on the evolution of the lateral 
distribution. 
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FIG. 12: Mean free path in air for the different muonic interactions as a function of the initial 
kinetic energy. This figure is courtesy of Sergio Sciutto. 



One of the earliest parameterizations of the muon LDF in vertical showers was empirically 
derived by Greisen [166]. This LDF was inspired by the NKG parametrization, and is 
factorized into two terms, 

p,{r) = N,{t) f,{r) , (36) 
where Nfj_{t) gives the normalization as a function of depth t, 

is a structure function describing the lateral shape of the shower, and tq = 320 m is analogous 
to the Moliere radius. 

Later, Vernov and collaborators proposed a semi-analytical form of the structure func- 
tion [187], 

/»-(^)" exp(-^) . (38) 

with F = 0.4 and tq = 80 m. Similar approaches were also suggested by Hillas' group at the 
University of Leeds [188] and by the SUGAR Collaboration [189]. The slopes from Eqs. (37) 
and (38) are in very good agreement with each other at intermediate distance, but Eq. (38) 
predicts a distribution which is flatter close to the shower core and more strongly damped 
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at large distances. These LDFs have been used to fit experimental data. However, neither 
function reproduces the whole radial range of an extensive air shower. This is a consequence 
of neglecting the shower age in formulating the structure functions. Very recently, the 
KASCADE Collaboration has used an NKG formula to fit muon density distributions [190]. 
The fits were performed close to the shower core (r < 200 m) with non-conventional values 
of vm and age parameter, s. 

One expects there to be a dependence of the LDF parameters on the shower age. However, 
in contrast to electrons, muons in an air shower are less attenuated and httle affected by 
Coulomb scattering, so the dependence of the LDF on the shower age is not the same as that 
exhibited by the electromagnetic component. The lateral growth of the shower is largely 
determined by the direction of emission of the parent particle and hence increases while 
the shower propagates downward. Two approaches for including shower age- dependence 
in the muon structure function have been discussed in the literature [167, 191]. Here, we 
consider the more recent treatment, in which a Vernov-like approach is used taking a slope 
dependence on atmospheric depth, F = 2 — s, with s as given in Eq. (34). 

It is easily seen [192, 193] that, if the parent particles are created with a p^, distribution, 
Pt/Po exp(— Pj,/po)'^Pt/Po, then the Vernov distribution at ground level has a value of ro 
given by 

ro = l iK) (39) 

where (p^) = 2pQ is the mean transverse momentum, (E^) the mean energy of muons and 
(hp) the mean height of muon production. These approximate expressions can serve to 
calculate the variation with depth of the parameters characterizing the lateral spread. The 
ratio (pj) I (Efi) can be considered constant while the shower develops [194] and the variation 
of To with altitude is determined only by the dependence of {h.p) on depth, t. 

Muons arc produced in every pion generation and their energy distribution follows that of 
their parents. There are three phenomena contributing to the behavior of (hp) as a function 
of t. The first is simply the dependence of the atmospheric density on height and tempera- 
ture. For an isothermal atmosphere of scale ho, one obtains (hp) oc ho \n{t/tp). The second 
phenomenon is the "pionization" process: the competition between pion production and 
decay. The last contribution to the behavior of (hp) is associated with systematics induced 
by hadronic interaction models. In what follows, we leave aside the issue of systematic errors 
and as an example adopt QGSJET 98 as the hadronic interaction model. Combining all these 
considerations, the characteristic radius ro{t) becomes. 



2 (p^) tGL , 



(40) 



where the subscript GL indicates that a quantity is given at ground level. For t — ^gl, 
Eq. (39) is recovered. 

In Fig. 10 we show the fi~^fi~ density distributions for a single lO^'' GeV proton shower at 
various depths. Fits to the Vernov-like distribution are overlaid on the simulation results, 
indicating validity of the parametrization [167]. Furthermore, the total number of muons 
from the fits agrees quite well with the corresponding values predicted by the Monte 
Carlo simulation, even though the fits are performed at core distances r > 100 m. 

Muons can travel long distances without interacting with the medium, and consequently 
the ground density profiles are significantly modified by the Earth's magnetic field. The 
global shower observables, hke longitudinal and lateral distributions, are not affected by the 
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geomagnetic field for zenitli angles d < 70° [181]. However, for the case of very inclined 
showers, which are dominated by muons, the density at ground is rendered quite asym- 
metric by the geomagnetic field. In the remainder of this section we describe these effects 
quantitatively [194]. 

Consider a highly relativistic muon of energy E^k. cp and transverse momentum that 
travels a distance d to reach ground. This muon suffers a deviation r from the shower axis 
given by 

r c ^ . (41) 

Now, it is easily seen that if the energy spectrum of muons is taken as 0(-E') = A E~'^ , the 
muon density is given by [194] 

P,{r) = ^{cp,df-^ T-^+^ . (42) 

To take into account the effect of the geomagnetic field, define a cartesian coordinate, 
{x,y), in the plane transverse to the shower axis, with y aligned to B±. The circular 
symmetry of the shower is distorted depending on B±, the distance traveled by the muons, 
and their energy distribution. For very large zenith angles the pattern results in two lobes 
on each side of the shower axis, one for the negatively and one for the positively charged 
muon components. The magnetic deviation 5x experienced by muons of different charges 
is [194] 

Sx = , (43) 

where e is the electron charge. Combining this with Eq. (41) we obtain, 

(^)0(<gv)"'— • <«) 

where f corresponds to the muon deviation in the transverse plane in the absence of a 
magnetic field, and a measures the ratio of displacement in the transverse plane due to 
as well as the displacement due to the magnetic field. The density of muons in the transverse 
plane can be obtained by making the transformation 



(46) 



X = X + ay + , y = y , (45) 

where the barred and unbarred coordinates indicate the position of the muon in the trans- 
verse plane in the absence and presence of the geomagnetic field, respectively. The muon 
number density reads 

oixv) 

P^.{x,y)=p^{x,y) 1^^^ 

where p^{x,y) is the density at a distance r = (x^ -|- y^Y^"^ in the case B — and the last 
factor is the Jacobian of the transformation. 

In a realistic situation the transverse position of the muon f is affected by both multiple 
scattering and the transverse position of the parent pions. Following [194], to account for 
this effect we use Eq.(41), setting d to a constant and we assume that at a given f there is 
an energy distribution. For convenience one introduces the variable e = log^g Efj_, such that 

(e) = A-7logior. (47) 
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FIG. 13: Lateral distribution of muons from AIRES simulations superimposed over the best fits 
obtained using Eq. (51). Prom top to bottom the curves correspond to 60°, 70° and 80° [194]. 



The muon density is taken to be 

p^(f,e) = P(e; {e),a) p^{f) , 



(48) 



where P is a distribution of mean (e) and standard deviation a. Now, one obtains the muon 
number density in the coordinate system {x,y), by using Eq. (46), 



where 



X — 



2E,, 



1/2 



(49) 



(50) 



The muon number density given in Eq. (49) depends on 3 quantities: (i) the distribution 
of e that hereafter is taken as a Gaussian with mean given by Eq. (47) and a ~ 0.4, (ii) 
the effective distance to the production point d, and (in) the lateral distribution function of 
the muons in the transverse plane. Figure 13 shows fits [194] to the lateral distributions at 
different zenith angles using the NKG-like LDF, 



p(r) — Mr 



1 + - 

n 



(51) 



with N, ip, K, TZ as given in Table I, and d taken as 16 km, 32 km, and 88 km, for 9 = 60°, 
70°, and 80°, respectively. One can see from the figure that the parametrization reproduce 
the simulation quite well up to a core distance of 1 km. 
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TABLE I: Best values for the parameters in Eqs. (47) and (51) as obtained from fits to Monte 
Carlo simulations, using SIBYLL 1.6 to process the hadronic interactions. 



9 


A 


7 


M 


V' 


K 


7^ 


60° 


2.67 ±0.23 


0.75 ±0.1 


569.9 


0.52 


4.05 


782.8 m 


70° 


4.04 ± 0.20 


0.73 ± 0.06 


227.1 


0.49 


4.35 


1010.0 m 


80° 


3.63 ±0.20 


0.81 ± 0.07 


78.4 


0.52 


4.49 


1513.0 m 



2. Muon content of the shower tail 

As discussed in the previous section, once the shower particle energies fall below eo, 
ionization losses take over from other electromagnetic processes, and the number of electrons 
and photons in the shower begins to decrease, while the number of muons remains more- 
or-less the same. We will refer to this region of the shower as the "tail." Most ground 
arrays are located below the altitude at which Xmax occurs, even in the case of vertical 
showers induced by ultra high energy primaries. This means that ground arrays observe 
predominantly shower tails. In this section we describe the variation of the shower tail's 
muon content with energy, and compare the original calculations of Hillas from the early 
1970's with more recent detailed Monte Carlo simulations. 

Muons are produced when a shower has cooled sufficiently to allow pions to decay before 
they interact (recall that the probability for decay of a 1 TeV pion in the atmosphere is 
somewhat less than 10%). Hillas' group at the University of Leeds reviewed various models 
for this cooling process and analyzed their consistency with data from emulsion experiments 
as well as cosmic air shower observations [195, 196]. They found the data at ground level 
to be best reproduced by the model (so-called "E") which predicted that for 9 = 14°, the 
number of muons in proton showers scales as 

oc . (52) 

Interestingly, this result differs only by an offset in the normalization when compared to the 
prediction from full-blown modern-day Monte Carlo simulations, as shown in Fig. 14. Of 
course, the exponent in Eq. (52) varies with zenith angle. It is possible to take into account 
the zenith angle dependence either through the "constant intensity cut" method [198], or 
by simply determining the behavior of the exponent as a function of zenith angle [182]. 
Furthermore, it has been recently noted that this exponent is more accurately taken to have 
a logarithmic energy dependence [86]. 

The muon content of EAS at ground level iV^, as well as the ratio N^/Ne, are sensitive 
to primary composition (here, is the electron content at ground level). To estimate 
the ratio of the muon content of nucleus induced to proton induced showers, we can resort 
again to the principle of superposition. Using Eq. (52) we find that the total number of 
muons produced by the superposition of A individual proton showers is, A^"^ oc A{E^/ A)^'^^. 
Consequently, in a vertical shower, one expects a cosmic ray nucleus to produce about 
more muons than a proton. This implies that a shower initiated by an iron nucleus produces 
about 27% more muons than a proton shower. Note, however, that a change in the hadronic 
interaction model could produce a much larger effect than a change in the primary species. 
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FIG. 14: Left panel: Total number of muons at ground level as a function of the shower energy. 
The dashed line indicates the Hillas parametrization for model "E", with a threshold energy set 
to that of the SUGAR experiment [26]. (For vertical showers considering the SUGAR'S energy 
threshold, one obtains an exponent 0.93 rather than the 0.94 used in the text.) The stars illustrate 
the results obtained from simulations carried out with aires + QGSJET 01, assuming proton 
primaries. The particles were injected vertically and the observation level was placed at 250 m 
above sea level. Muons with energies below the threshold 0.75 GeV are not taken into account in 
the simulations [183]. If the hadronic interactions are modeled with sibyll 1.6 rather than QGSJET 
01, an exponent of 0.88 best fits the simulation [197]. Right panel: Ratio of the muon content for 
EAS produced by primary gammas and protons. The geomagnetic field is set to the PAO Southern 
site [150]. 



For example, replacing QGSJET 01 with sibyll 1.6 as the hadronic interaction model leads 
to a prediction of 60% more muons for an iron shower than for a proton shower [182]. 

The situation for gamma-induced showers is a bit different. In this case the muon com- 
ponent of the shower does not simply follow Eq. (52) because of the LPM and geomagnetic 
field effects [150]. Competition between the two processes leads to a complex behavior in 
N'^/NJ^, as shown in Fig. 14. 

In this section we have described the four main energy loss mechanisms for muons en route 
through the atmosphere. The rate of energy attenuation for muons is much smaller than it is 
for electrons, and the energy loss processes are only really of interest in the case of extremely 
inclined showers for which the original electromagnetic component is mostly absorbed. In 
such cases, small electromagnetic sub-showers can still arise from bremsstrahlung, pair pro- 
duction and knock-on electrons. In addition, muon-nucleus interactions induce hadronic 
sub-showers. We also discussed the effect of different energy loss mechanisms on the elec- 
tron and muon distributions in time and space. Because they are less subject to multiple 
scattering, muons tend to arrive at the ground earlier and more compressed in time. The 
ratio of muons to electrons far from the core is much greater than it is near the core, and 
this effect is more pronounced at higher zenith angles. 
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The muon content of the shower tail is quite sensitive to unknown details of hadronic 
physics. This implies that attempts to extract composition information from measurements 
of muon content at ground level tend to be systematics dominated. The muon LDF is 
mostly determined by the distribution in phase space of the parent pions. However, the 
pionization process together with muon deflection in the geomagnetic field obscures the 
distribution of the first generation of pions. A combination of detailed simulations, high 
statistics measurements of the muon LDF and identification of the primary species using 
uncorrelated observables could shed hght on hadronic interaction models. 

IV. FINGERPRINTS OF PRIMARY SPECIES IN EAS 

A determination of primary composition is invaluable in revealing the origin of cosmic rays 
as this information would provide important bounds on sources and on possible production 
and acceleration mechanisms. In addition, a proper interpretation of anisotropy information 
requires knowledge of the primary mass due to the influence on propagation of the galactic 
and intergalactic magnetic flelds. Attempting to determine the primary composition of 
cosmic rays by measuring various shower parameters is fraught with systematic uncertainty. 
Furthermore, because of the stochastic nature of the extensive air showers, there are inherent 
shower-to-shower fluctuations in measured shower observables that cannot be attributed to 
experimental systematic error alone. Therefore, the determination of primary composition 
on an event-by-event basis is an intractable problem. Nevertheless, statistical analyses of 
shower observables known to correlate with the primary composition are possible. Based on 
the general signatures of the EAS described in previous sections, we provide a summary of 
the observables that help to separate primary species. 

A. Photon showers 

In this section we provide an overview of how the EAS characteristics described in the 
previous sections allow one to distinguish photon primaries from other species. As discussed 
in section III A photon-induced showers are expected to generate fewer muons than baryon- 
induced showers. This clear signature can be exploited by surface arrays which are equipped 
with dedicated muon detectors or arc capable of distinguishing muons using shower observ- 
ables sensitive to the muon content. The AGASA Collaboration [199, 200] has used the 
muon content of the detected EAS to set bounds on the percentage of photon primaries 
present in the observed flux. AGASA comprises 111 stations covering an area of 100 km^. 
Detectors of 2.8 — 10 m^ area, capable of measuring muon densities up to ^ 10 m~^, were 
deployed in 27 stations. The analysis of the AGASA Collaboration, which takes into ac- 
count the LPM effect and conversion in the geomagnetic field, shows that at the 95% CL, 
the fraction of 7-rays above 10^° GeV, lO^^-^s GeV, and lO^^-^ GeV is less than 34%, 59%, 
and 63%, respectively. Of course, these bounds depend on the hadronic interaction models 
used to simulate the showers. Several models were used in this analysis, and the reported 
limits are the least restrictive ones. 

Another powerful tool for discriminating between photons and baryons using data col- 
lected by surface detectors relies on comparing the flux of vertical showers to that of inclined 
showers, a technique which exploits the attenuation of the electromagnetic shower compo- 
nent for large slant depths. As an illustration of this technique, we describe below the 
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FIG. 15: Integral number of inclined events as a function of energy for the Haverah Park data 
set compared with the predictions for iron, proton, and photon primaries. Here the energy is 
calculated assuming a proton primary. The slope of the assumed primary spectrum oc E~^-'^^ is 
shown to illustrate the increase of triggering efficiency with energy [201]. 

constraints on the gamma ray flux obtained from Haverah Park measurements [201]. The 
first crucial ingredient in the analysis is the vertical flux normalization. This should be 
determined in a way which is free from systematic uncertainties associated with the primary 
composition. Fluorescence detectors, which record "calorimetric" measurements, provide 
the best tool to attain this normalization, and in the analysis described here the data from 
Fly's Eye [41] were used [202]. From this known vertical spectrum and Monte Carlo sim- 
ulations of shower propagation and detector response, a prediction can be made for the 
expected rate of inclined events for each type of primary. For inclined showers in the zenith 
angle range 60° < 9 < 80°, the Haverah Park experiment collected 46 events with energy 
above 10^" GcV and 7 events above 10^°'^ GeV. A Comparison of these observations to the 
results extracted from simulations is shown Fig. 15. If one assumes the primary spectrum 
comprises a mixture of protons and photons, then the Haverah Park data imply that above 
10^° GeV, less than 30% of this admixture can be photons, and above 10^°-^ GeV less than 
55% can be photons. Both of these statements are made at the 95% GL [201]. Even though 
Fly's Eye provides a flux measurement which is independent of the mass composition, one 
should be aware of the inherent systematic uncertainties in aperture estimates of fluores- 
cence detectors. A separate normalization technique using both fluorescence and surface 
array data leads to bounds within 20% of the previous estimates [197]. 

The sensitivity of PAG for isolating gamma ray primaries using this method was estimated 
in [203] . Given the huge statistics - above 60° the aperture of the observatory is increased 
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FIG. 16: Distribution of observed X^ax from HiRes stereo data for showers in the energy range 
10^'^ — 10^'^ '^ GeV (solid histogram). The predictions for proton (upper figure) and for iron (lower 
figure) are given for QGSJET 01 (dashed histogram) and sibyll 2.1 (dotted histogram) This figure 
is courtesy of Greg Archbold. 



by almost a factor 2 - PAO will place severe constraints on the photon content of the 
observed flux. Additionally, hybrid techniques available to PAO will facilitate independent 
cross-checks on the systematic uncertainties. 

The effect of the geomagnetic field on photons also leads to an energy dependence and 
characteristic anisotropy in the directional distribution of the fraction of events with abnor- 
mally deep profiles which are not compatible with proton or nuclei primaries. This technique 
has not yet been implemented in the analysis of real data, but the potential for the HiRes 
and PAO experiments has been evaluated [154, 155]. 

B. Hadronic primaries 

We now discuss how baryonic species may, to some extent, be distinguished by the signa- 
tures they produce in the atmosphere. As in the previous section, we consider both surface 
array and fluorescence detector observables. 

As mentioned in Sec. Ill A, the main purpose of fluorescence detectors is to measure 
the properties of the longitudinal development. The shower longitudinal profile can be 
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parameterized using the Gaisser-Hillas function [204] 



V V \ [(^max~-^l)/A] f V V 

^ — ^1 \ I ^max — ^ 



7Ve(X) = 7Ve,^ax (^^— 3^ j ^""P \ /' ^^^1' (^^j 

where iVe,max is the size at the maximum, Xi is the depth of the first observed interaction, 
and A is a floating parameter in the fit, generally fixed to 70 g/cm^. Using this parametriza- 
tion, fluorescence detectors can measure X^^x with a statistical precision typically around 

30 g/cm^. 

The speed of shower development is the clearest indicator of the primary composition. It 
was shown in Sec. Ill A using the superposition model that there is a difference between the 
depth of maximum in proton and iron induced showers. In fact, nucleus-induced showers 
develop faster, having Xmax higher in the atmosphere. From Monte Carlo simulations, one 
finds that the difference between the average Xmax for protons and iron nuclei is about 90 - 
100 g/cm^. However, because of shower-to-shower fluctuations, it is not possible to obtain 
meaningful composition estimates from Xmax on a shower-by-shower basis, though one can 
derive composition information from the magnitude of the fluctuations themselves. For 
protons, the depth of flrst interaction fluctuates more than it does for iron, and consequently 
the fluctuations of X^ax sue larger for protons as well. Speciflcally, the standard deviation 
c(-^max) is 53 g/cm^ for protons and 22 g/cm^ for iron [9]. These fluctuations depend only 
weakly on the choice of interaction model. The HiRes Collaboration has recently analyzed 
their stereo data sample in the energy range E = 10^'^ — lO^*^'^ GeV [205]. The results are 
shown in Fig. 16, together with the expected distributions using several hadronic interaction 
models. While agreement between data and Monte Carlo is not perfect for any of the models, 
the data do qualitatively suggest a proton dominated composition. 

Changes in the mean mass composition of the cosmic ray flux as a function of energy 
will manifest as changes in the mean values of Xmax- This change of Xmax with energy is 
commonly known as the elongation rate theorem [206]: 

For purely electromagnetic showers, Xmax(-E) ~ Xq ln{E/eo) and then the elongation rate 
is DgRi Xq. For proton primaries, the multiplicity rises with energy, and thus the resulting 
elongation rate becomes smaller. This can be understood by noting that, on average, the 
flrst interaction is determined by the proton mean free path in the atmosphere, A at. In this 
flrst interaction the incoming proton splits into {n{E)) secondary particles, each carrying 
an average energy E/(n{E)). Assuming that X-^^{E) depends logarithmically on energy, 
as we found with the Heitler model described in Sec. Ill A, it follows that, 

Xr^UE) = Ajv + Xo ln[^/(n(^))] . (55) 

If we assume a multiplicity dependence {n{E)) w uqE^, then the elongation rate becomes, 

5\n{n{E)) 



5\nE 



which corresponds to the form given by Linsley and Watson [207] , 



De^Xo 



6ln{n{E)) A7v51n(Ajv) 
6\nE 5\nE 



Xo {1-B) . (57) 
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FIG. 17: Left panel: Variation of Xmax with energy as seen by different experiments: Fly's 
Eye [208], HiRes-MIA [209, 210], HEGRA-AIROBICC [211], CASA-BLANCA [212], DICE [213], 
SPASE- VULCAN [214], and YAKUTSK [215]. The rising curves indicate simulated results for 
proton and iron primaries using various hadronic interaction models [162]. Right panel: The cir- 
cles indicate the experimental measurements of log at 1000 m from the core vs. logarithm of 
the primary energy. The lines indicate the la band for iron, proton, and photon predictions [199]. 



Using the superposition model introduced in Sec. Ill A and assuming that 

is not changing with energy, one obtains for mixed primary composition [207] 

9(lnA)' 



= Xo (1 - B) 



dlnE 



(59) 



Thus, the elongation rate provides a measurement of the change of the mean logarithmic 
mass with energy. One caveat of the procedure discussed above is that Eq. (57) accounts 
for the energy dependence of the cross section and violation of Feynman scaling only for 
the first interaction. Note that subsequent interactions are assumed to be characterized 
by Feynman scaling and constant interaction cross sections, see Eq. (58). Above 10'' GeV, 
these secondary interactions play a more important role, and thus the predictions of Eq. (59), 
depending on the hadronic interaction model assumed, may vary by up to 20% [86]. 

For convenience, the elongation rate is often written in terms of energy decades, Diq = 
d{Xjaax) /dlogE, where Diq = 2.30^. In Fig. 17 we show the variation of (X^ax) with 
primary energy as measured by several experiments together with predictions from a variety 
of hadronic interaction models. For protons and iron, Monte Carlo simulations indicate 
DiQ ?a 55 g/cm^ and for photons Diq ^ 84 g/cm^ [9]. Recent results presented by the 
HiRes Collaboration [205] using stereo data favours a light component above 10^ GeV, 
and the reported variation of (Xmax) with logarithm of primary energy is Diq = 54.5 ± 6.5, 
consistent with a constant or slowly changing composition between 10^ GeV and lO^"'^ GeV. 
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FIG. 18: Comparison of the simulated distribution of r] (histogram) with measured distributions 
(points) for iron (left) and proton (middle). One can see that an iron-dominated composition best 
fits the data, but that the addition of a lighter component is needed to fit the distribution at 
large values of rj. The right panel shows the best composition fit to the measured distribution of 77 
(points) from maximum likelihood analysis. 



As an attentive reader should know by now, the muon content of the shower is strongly 
sensitive to the nature of the primary. The AGASA Collaboration used measurements of 
the muon component to discern the primary composition [199]. For events with estimated 
energy > 10^° GeV, and zenith < 36°, the muon density at 1000 m from the shower core was 
used to estimate the primary mass. Expected muon densities for iron and proton primaries 
were estimated from Monte Carlo simulations, and comparison of the expected to observed 
densities suggests a dominance of light composition. Specifically, above 10^*^ GeV the average 
fraction of iron is 14^^^%, rising to 30l^% above 10^°-^^ GeV, and above 10^°-^ GeV the 
fraction is found to be below 66% at the la level. In Fig. 17 we show the distributions of 
muon density at 1000 m from the core as a function of primary energy as reported by the 
AGASA Collaboration, together with the predictions for ±lcr bounds for iron, proton and 
photon induced showers. 

The steepness of the lateral distribution of particles at ground level also correlates with 
the depth of shower maximum, and thus carries information about the primary species. For 
instance, a proton-initiated shower would have a steeper average lateral distribution, since 
the shower develops deeper in the atmosphere than an iron-initiated shower of the same 
energy. Recently, this approach has been used in conjunction with the latest shower and 
detector simulation codes to re-interpret the lateral distribution measurements from Hav- 
erah Park [216] and Volcano Ranch [217-219]. In the case of the Volcano Ranch array, 80 
scintillators were laid out in a grid with a separation of 147 m, facilitating a very fine-grained 
measurement of the lateral distribution. Recall that the signal measured by plastic scintil- 
lators is the average energy loss in the scintillator of electrons, muons and photons in units 
of the energy loss of vertically penetrating muons. The correpsonding lateral distribution 
was parameterized by an NKG-like formula [168], 

where N is the shower size, rj and a describe the logarithmic slope, and tm — 100 m at 
Volcano Ranch. For events at median energy 10^ GeV and shower sizes A'" = 4xl0^ — 6x10^, 
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FIG. 19: The triangles (iron) and circles (proton) indicate the rise time (left) and fall time (right) 
as a function of C- Solid symbols correspond to a primary zenith angle of 35°, while the open 
symbols correspond to 60° [227]. 

measurements of r] (with a fixed to 1) were analyzed in the range of zenith angle 1.0 < 
sec^ < 1.1 [220, 221]. These measurements are shown in Fig. 18 along with the recently 
simulated results for purely proton and iron primaries using qgsjet 98 as the hadronic 
interaction model. One can clearly see the dependence on primary composition reflected 
in the distribution of rj. To quantify this dependence, a maximum likelihood technique 
was employed assuming a bimodal composition of proton and iron. The cosmic ray mass 
composition, deduced from Volcano Ranch data, resulted in a mean fraction (89 ± 5(stat) ± 
12(sys)) % of iron in the whole energy range 10^''^ GeV to lO^'' GeV; mean energy 10^ GeV. 
The resulting admixture is also shown in Fig. 18. Systematic uncertainties associated with 
the hadronic interaction model are important in this analysis. If QGSJET 98 is replaced by 
QGSJET 01, then the fraction of iron changes from (89 ±5)% to (75 ± 5)%. 

As described in Sec. Ill A the electromagnetic component of an EAS suffers more scat- 
tering and energy loss than the muonic component and consequently, muons tend to arrive 
earlier and over a shorter period of time. This means that parameters characterizing the 
time structure of the EAS will be correlated with X^ax and hence with primary mass. An 
early study of the shower signal observed in water Cerenkov detectors at the Haverah Park 
array [222] established the utility of a shower property known as rise time in estimating the 
primary composition. Specifically, the rise time, ti/2, is defined as the time for the signal 
to rise from 10% to 50% of the full signal. Interestingly, the relation between rise time and 
depth of maximum also implies a relation between the rise time and the elongation rate. As 
suggested by Linsley [206] if P represents the average value of some shower parameter, such 
as the rise time, which does not depend explicitly on primary energy but depends on the 
depth of observation, X, and the depth of shower maximum Xmax, then the elongation rate 
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FIG. 20: Time distributions from data collected by the PAO during the period of May to November 
2002 in the radial range 600 — 1400 m [227]. In each plot the upper points corresponds to 1 < 
sec 6 < 1.3 while the lower corresponds to 1.3 < sec 9 < 2. A fit to a linear cosine is overlayed on 
the points, where the fitted parameter "off" is the mean value and "amp" is the amplitude of the 
asymmetry [227]. 



can be derived from the following expression: 

^). — 

where F depends on the depth dependence of P. For a depth dependence of the form 

f{X/Xmax), F = X/Xmax, wheieaS for f{X - Xmax), F =1. 

This alternative approach was applied to Haverah Park data using the rise time ti/2 of the 
signals [223] to produce a measurement of the elongation rate from 10^'^ GeV to 10^^ GeV. 
By means of an experimentally determined value for the dependence of ti/2 with depth, they 
obtained an elongation rate of 70 ±5 g/cm^, averaged over the previously mentioned energy 
range. Their result is suggestive of an evolution to lighter species with rising energy. 

Recently, an extension of the work on the shower front thickness using Haverah Park data 
was performed, focusing on the highest energy events [224]. In this analysis the averaged 
values of the rise time at a large distance from the core were compared with Monte Carlo 
(corsika/qgsjet 01) predictions for proton and iron values. The result indicates a large 
fraction 80%) of iron nuclei at 10^° GeV [225]. 

Azimuthal asymmetries in the size [226] and time structure of signals at the ground [227] 
have been observed in non-vertical showers. The origin of these asymmetries has been 
discussed in Sec. Ill A. The AUGER Collaboration has found that the asymmetry in time 
distributions offers a new possibility for the determination of mass composition, because its 
magnitude is strongly dependent on the muon to electromagnetic ratio at the observation 
level. A preliminary study of the timing information of EAS using simulated proton and 
iron events was used to estimate the sensitivity of the PAO in discrimination of baryonic 
primaries [227]. The following observables were analyzed: the rise time (time between 10% 
and 50% of the integrated signal), fall time (time between 50% and 90%) and the time 
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FIG. 21: Iron fraction from various experiments: Fly's Eye (A), AGASA AlOO (■), AGASA Al (□) 
using SIBYLL 1.6 ([228] and references therein) and Haverah Park [216], using QGSJET98 (o)). The 
mean composition determined in [217] with the corresponding error for the Volcano Ranch energy 
range using QGSJET98 (*) is shown. The solid line arrow indicates the recent result using rise time 
measurements from Haverah Park [224]. The dashed arrow lines represents upper limits obtained 
by the AGASA Collaboration with QGSJET98 [199]. The dot dashed horizontal line corresponds 
to results reported by the HiRes Collaboration [205] 



between 10% and 90% of the signal. The timing variables as a function of the azimuth angle 
in the shower plane, (, at fixed range of core distances for proton and iron showers are shown 
in Fig. 19. The incoming direction of the shower corresponds to C = 0. As one can see from 
the figure, these distributions tend to flatten with increasing primary mass. A first analysis 
seems to indicate that the fall time would be a better discriminating tool. One expects 
stronger asymmetries at intermediate core distances, where the electromagnetic component 
dominates. In Fig. 20 we show the mean rise time and fall time as a function of (, for events 
with energy above 10^ GeV collected by the PAO in the radial range 600 — 1400 m. These 
results, while preliminary, indicate the promise of this method for composition studies, once 
large statistics samples become available. 

In summary, the primary composition has been studied over various energy ranges using 
several experimental techniques. The variety of results is summarised in Fig. 21, which shows 
the fraction of iron as a function of energy. Surface arrays such as Haverah Park, Volcano 
Ranch and Akeno- AGASA infer X^ax; and hence the overall composition, from properties of 
secondary particle distributions at the ground. In this case, the primary source of systematic 
error arises from uncertainties in the hadronic interaction models. Fluorescence detectors 
such as Fly's Eye and HiRes observe an image of the longitudinal shower profile and extract 
-'^max directly. Consequently such measurements do not suffer from uncertainties in hadronic 
event generators, though the data analysis still faces the challenge of assessing atmospheric 
properties as a function of time. Future stereo data from HiRes and hybrid PAO observations 
will provide a higher statistics sample with a better control of the systematic uncertainties 
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and will certainly provide clues to the nature of ultra high energy cosmic rays. 

V. DEEPLY PENETRATING SHOWERS 

Even at large zenith angles, the mean free path of a neutrino in the atmosphere is much 
larger than the atmosphere's slant depth. However, nearly horizontal showers are especially 
interesting since in this case the likelihood of interaction is maximized, and furthermore, 
neutrino induced showers can be easily distinguished from those induced by hadrons high 
in the atmosphere. In this section we first consider strategics for detecting these kinds of 
signatures. After that, we focus attention on recent novel scenarios with TeV-scale quantum 
gravity and discuss the observables of neutrino showers associated with both sub-planckian 
and trans-planckian physics. 

A. Everyday weakly interacting neutrinos 

Neutrinos are unique and thus far relatively untapped astronomical messengers [229]. Up 
to now, the only directly observed extraterrestrial neutrinos are those from the Sun [230] and 
from supernova SN1987A [231, 232]; these are low energy (MeV-range) neutrinos. Higher 
energy neutrinos should be generated by cosmic "beam dumps" in which baryonic particles 
collide with interstellar media. These neutrinos are particularly appealing for astronomy 
since they are undefiected by magnetic fields and they can travel cosmological distances 
without interacting [233, 234]. In addition to providing a new window for astronomy, cosmic 
neutrino observations may also enlighten our understanding of fundamental physics. For 
instance, it will be possible to test Lorentz invariance at extremely high energies and to 
hunt for exotic processes such as neutrino decay, CPT violation, and small 5im? oscillations 
into sterile neutrinos (see e.g., [235-241]). Since neutrinos interact only weakly, very large 
detector volumes are required, ideally on the order of Ikm^we [242]. As discussed in the 
introduction the PAO overlooks ISkm^we [33], which is sufficient to collect a statistically 
significant sample of neutrino showers, provided they can be separated from the hadron and 
photon-induced showers. In what follows we discuss the characteristics of neutrino-induced 
EAS, and comment on the qualities which may provide a handle to separate these showers 
from background. 

In a typical collision in the Earth's atmosphere the neutrino (with energy Eij) scatters off a 
proton either via the charged current, {i'i,Vi)N — >• {l~ , -|- anything, or the neutral current, 
{i'i,i7i)N — > {vi,Vi) + anything. ^° The kinematics of these reactions can be characterized by 
the inelasticity parameter y — (1 — cosi?*)/2 and the 4 momentum fraction of the proton 
carried by the struck quark x = Q'^/{ys), where —Q^ is the invariant momentum transfer 
between the incident neutrino and the outgoing Icpton. For a given the lowest x is 
achieved when y = I and the lowest y when x = 1. Thus, kinematically the small values of 
X are associated with large values of y and vice versa. 



The rate of interaction of fg, i^fi, Vt, I'n, T^t, with atmospheric electrons is neghgible compared to inter- 
actions with nucleons. The case of interactions is exceptional because of the intermediate Glashow 
resonance formed via W boson production at -E^^ « 10^'^ GeV [243]. 
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The charged current differential cross section of a neutrino scattering on an isoscalar 
nucleon N = {n + p)/2 of mass M is given by 



2^CC 



uN 



dx dy 



TT 



1 + (1 - yf 



(62) 



where mw ^ 80.423 GeV denotes the W boson mass, Gp = 1.16639 x 10"^ GeV^ is the 
Fermi coupUng constant and the structure functions Fj in terms of the quark distribution 
functions of the nucleon qi{x, Q^) read 



F^^ = 2x 



+ ds + Us + Sg + Cs + bs + t, 



xF^'^ = 2x {ds + Ss-Us- Cs) , 



(63) 
(64) 



and 



^cc ^ ^cc _ 2^fCC 



L --2 ---1 ■ (65) 

Here the subscripts v and s label valence and sea contributions, and m, d, c, s, b denote 
the distributions for various quark flavors in a proton. In the deep inelastic factorization 
scheme, Eq. (62) can be re-written in terms of quark distributions as [244, 245] 



d^a^^2G^ME^ 
dxdy TT 



m 



w 



^ [xq''''M') + {l-yfxt''M^)\ , (66) 



where 



CC/ /-)2\ '^v{,^i Q^) + dy{x^ Q'^) Us{x, Q"^) + ds{x, Q'^) ^ ^ ^ ^ , i ^ ^2^ 



+ Ss{x,Q^)^bs{xM^) , (67) 



and 



+ c,(x,(5^) +t,(a;,Q^) 



(68) 



In Eq. (66) we have omitted perturbative QCD corrections, which for \/s > 10^ *^ GeV are 
insignificant. The average energy loss of this process is (y) ~ 0.2. Duplicating this procedure, 
one can straightforwardly obtain the cross section for the neutral current neutrino-nucleon 
interaction. 
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dx dy 27r + 

where the quark densities are given by 

'uyix.Q"^) + dy{x,Q'^ 



[xq''^{x,Q') + {l-yrxt''{x,Q')] , 



(69) 



+ 2 



Wvf + {9\f + {9lf + {9\f\ 



+ 2[.,(a;,g2) + 6.(a;,g2)] [(<?^)^ + (<?!)] 
+ 2[c,(x,Q^) + i,(x,g^)][(^^)^+(^l)2] , 



(70) 
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Us(x, Q"^) + ds{x, g^ 
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+ 2[.,(x,g2) + 6,(x,g2)] + (4)] 
+ 2[c,(x,g2) + i,(x,g2)][(^;>r + (^^f] 



(71) 



Here, mz ~ 91.187 GeV is the mass of the neutral intermediate boson and 

1 2 . o„ .1 

2 + 3' 



9a- 2' 



(72) 
(73) 



4 

2 - sin^ , " 2 ' 

are the vector and axial-vector couplings for down- and up-type quarks, respectively; with 
sin^^w — 0.226 the weak mixing parameter [137]. Similar calculations lead to the cross 
sections for VN scattering. For further details see e.g., [246]. 
The X — region probed by ultra high energy neutrinos. 
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(74) 



is well beyond the region accesible by the HERA experiments (see Fig. 3). As we discussed in 
Sec. II, in the renormalization group-improved parton model, the structure functions are ex- 
trapolated to very low x considering leading order (LO), next to leading order (NLO), and/or 
double-leading-logarithmic evolution of DGLAP equations. Using the CTEQ4 pdf 's [247] 
one finds [248]: 



and 
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(75) 
(76) 
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For lO'' GeV ^ E^ < 10^^ GeV, these cross sections are correct to about 10%, which is 
smaller than the systematic uncertainties that cosmic ray experiments typically contend 
with. Note that the reason this uncertainty is small compared to the uncertainty in the 
cross section shown in Fig. 5 is a consequence of the existence of two viable models for cross 
section extrapolation in the case of pp interactions. Neutrino interaction lengths 



1.7 X 10^ kmwe 



pb 



(CN)C 



(79) 
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are therefore far larger than the Earth's atmospheric depth, which is only 0.36 kmwe even 
when traversed horizontally. As a consequence, neutrino showers, unlike baryon or photon 
induced showers, can begin deep in the atmosphere. So, to obtain a clean signal of neutrino- 
induced showers one should be able to identify deeply developing cascades in the whole 
sample of EAS. 

For large zenith angles {6 > 70°), an air shower initiated by a neutrino can be distin- 
guished from that of an ordinary hadron by its shape. As discussed in Sec. Ill, baryons 
interact high in the atmosphere. Consequently, at ground level the electromagnetic part of 
these inclined showers is totally extinguished (more than 6 equivalent vertical atmospheres 
were gone through) and only the muon channel survives. Besides, the shower front is ex- 
tremely flat (radius of curvature > 100 km) and the particle time spread is very narrow 
{At < 50 ns). Since neutrinos can interact deeply in the atmosphere, they can initiate 
showers in the volume of air immediately above the detector. Therefore, quasi-horizontal 
showers initiated by neutrinos would still present a curved front (radius of curvature of a 
few km), with particles well spread over time, (9(//s), allowing a good characterization of 
the cascade against background. 

If the primaries are mainly electronic and muonic neutrinos, as expected from pion decays, 
two types of neutrino showers can be distinguished: "mixed" (with full energy) or "pure 
hadronic" (with reduced energy), respectively [249]. In the charged current interaction of 
a z/g, an ultra high energy electron having about 80% of the z/g energy is produced and 
initiates a large electromagnetic cascade parallel to the hadronic cascade. In contrast, the 
charged current interaction of a i^^ produces a muon which is not easily detectable by existing 
experiments. In the presence of maximal z/^/z/^-mixing, i/T--showers must also be considered. 
The r mean flight distance is ~ 50E^/{10^ GcV) km, whereas the distance between the 
position of first impact and ground is ~ 30 km, hence only r's with energy < lO*'^ GeV will 
decay before reaching the ground. Since the r is produced with about 80% of the Ur energy, 
showers initiated by ultra high energy u^s will be indistinguishable from showers. 

Another interesting category of neutrino-related showers comprises events in which a 
neutrino skims the Earth, traveling at a low angle along a chord with length of order its 
interaction length. Some of these Earth-skimmers may be converted into charged leptons 
in the Earth's crust. Unhke electrons that do not escape from rocks, at the energies of 
interest, muons and tau leptons travel up to C(10 km) inside the Earth. Although muons 
do not produce any visible signal in the atmosphere, taus can produce clear signals for 
both fluorescence eyes [250] and surface arrays [251] if they decay above the detector. The 
phenomenon would thus increase the i^-event rate and enhance the sensitivity of the PAO 
to neutrino fluxes. 

Up to now we have only considered signals one might expect from Standard Model (SM) 
processes. Many scenarii with new physics beyond the electroweak scale, M^w, have been 
proposed, some of which increase the weak interaction cross section (see e.g., [252-255]), 
and hence would have observable implications. As an example of such SM extensions, we 
consider in the last part of this review a scenario which has generated a great deal of recent 
interest. 

B. Neutrino interactions mediated by gravity 

A promising route towards reconciling the apparent mismatch of the fundamental scales of 
particle physics and gravity is to modify the short distance behavior of gravity at scales much 
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larger than the Planck length. This can be accomplished in a straightforward manner [256- 
258] if one assumes that the SM fields are confined to a 4-dimensional world (corresponding 
to our apparent universe), while gravity lives in a higher dimensional space. One virtue 
of this assumption is that very large extra dimensions are allowed without conflicting with 
current experimental bounds, leading to a fundamental Planck mass much lower than its 
effective 4-dimensional value. In particular, if spacetime is taken as a direct product of 
a 4-dimensional spacetime and a flat spatial n-dimensional torus T"' (of common linear 
size 27rrc), one obtains a definite representation of this picture in which the effective 4- 
dimensional Planck scale, Mpi ~ 10^^ GeV, is related to the fundamental scale of gravity, 
Md, according to M|j = 87r Mj^^^r" [256]. Now, a straightforward calculation shows that, 
for n = 1, low-scale gravity within toroidal compactifications is excluded, as gravity would 
be modified at the scale of our solar system. Astrophysical constraints require Mjj ^ 10 TeV 
for n — 2,3 and Md > 4 TeV for n = 4 [259]. For n > 5, however, Md may be as low as a 
TeV [260-264]. 

From our 4-dimensional point of view the massless graviton appears as an infinite tower of 
Kaluza-Klein (KK) modes, of which the lowest is the massless graviton itself, but the others 
are massive. The mass squared of each KK graviton mode reads, — Yl^=i ^l/^l-: where the 
mode numbers are li e Z. Note that the weakness of the gravitational interaction is partially 
compensated by the tower of KK modes that are exchanged: the square coupling M^^ of the 
graviton vertex is exactly cancelled by the large multiplicity of KK excitations ~ s"/^ r", so 
that the final product is ~ s"/VM^+". Indeed, if one includes in the interaction the brane 
Goldstone modes, a form factor introduced at each graviton vertex [265] . This 

exponential suppression, which parametrizes the effects of a finite brane tension, provides 
a dynamical cutoff in the (otherwise divergent) sum over all KK contributions to a given 
scattering amplitude. Altogether, one may wonder whether the rapid growth of the cross 
section with energy in neutrino-nucleon reactions mediated by spin 2 particles carries with 
it observable deviations from SM predictions.^^ 

A simple Born approximation to the elastic neutrino-parton cross section (which underlies 
the total neutrino-proton cross section) leads, without modification, to (Td ~ [266, 267]. 
Unmodified, this behavior by itself eventually violates unitarity. This may be seen either 
by examining the partial waves of this amplitude, or by studying the high energy Regge 
behavior of an amplitude Ar^s, t) dtli spin-2 Regge pole, viz., intercept a{0) — 2. 

For the latter, the elastic cross section is given by 

dOel ^ IMsM ^ ^2a(0)-2 ^ ^2 /ggX 

di S2 ' ^ ^ 

whereas the total cross section reads as 

9m[Afl(s,0)] .a(0)-i - /eiN 
(7tot ' ^ ~ s, (81) 

s 

so that eventually, o"ci > o"tot [268]. Eikonal unitarization schemes modify these behaviors. 
Specifically, for large impact parameter, a single Regge pole exchange amplitude yields 
o^tot ~ ln^(s/so) [269, 270]. For short impact parameters, partial wave unitarity is a tall 



In what follows we only take into account KK excitations on the gravity sector without considering string 
recurrences of any other field. For a treatment of the latter see e.g., [254]. 
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order as corrections to the eikonal amplitude arc expected to become important. Note that 
graviton self interactions carry factors of t associated to the vertices, and thus as t increases, 
so too does the attraction among the scattered particles. Eventually it is expected that 
gravitational collapse to a black hole (BH) will take place, absorbing the initial state in such 
a way that short distance effects will be screend by the appearance of a horizon [271-274].^^ 
According to the Thome's hoop conjecture [276], a BH forms in a two-particle collision 
when and only when the impact parameter is smaller than the radius of a Schwarzschild 
BH of mass equal to V§ = ^/xs. The conjecture thus predicts a total cross section for BH 
production proportional to the area subtended by a "hoop" 

aBH = F{n)7rrl{^) (82) 

of radius [277, 278] 

"2"7r("-3)/2r(^ 



D 



M, 



D 



1 

1 + n 



n 



(83) 



where F{n) is a constant of order unity. Recent progress has confirmed the validity of 
Eq. (82) and evaluated the dimension-dependent constant F{n), analytically in four dimen- 
sions [279] and numerically in higher dimensions [280]. Of course, this work is purely in the 
framework of classical general relativity, and is expected to be valid for energies far above the 
Planck scale, for which curvature is small outside the horizon and strong quantum effects are 
hidden behind the horizon. Extending it to the planckian regime of center-of-mass energies 
close to M^) requires a better understanding of quantum gravity than we now possess. Thus 
it is important to impose a cutoff on the mass of microscopic BHs for which Eq. (82) can 
reasonably be expected to hold. 

In the course of collapse a certain amount of energy is radiated in gravitational waves by 
the multipole moments of the incoming shock waves, leaving only a fraction y = Mbh/V^ 
available to Hawking evaporation, where Mbh is a lower bound on the final mass of the 
BH. This ratio depends on the classical impact parameter b, and so the inclusive production 
of BHs proceeds through different final states for different values of b. These final states 
are characterized by the fraction y{z) of the initial parton center-of-mass energy which is 
trapped within the horizon. Here, z = b/b^ax: and biaax/fsi^^) is given in Table II [280]. 
With a lower cutoff MBH,min on the BH mass required for the vahdity of the semi-classical 
description, this implies a joint constraint 



> MBH,min (84) 



on the parameters x and z. Because of the monotonically decreasing nature of y{z), Eq. (84) 
sets an upper bound z{x) on the impact parameter for fixed x. The corresponding parton- 
parton BH cross section is d-^^{x) = nb'^{x), where b — zb„iax- The total BH production cross 
section is then [264] 

f^.iV^BH(-^i^>-^BH,min,MD) = / 2 dx fi{x, Q) a^^^ix) , (85) 

/ BH,min ' ^ 



This paper does not purport to be a comprehensive review of TeV-scale gravity BHs; for an up-to-date 
and detailed discussion of the topic, the reader is referred to [275] . 
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TABLE II: Gravitational collapse parameters. 



n 


2 


3 


4 


5 


6 


7 


VaxAs(\/^) 

F{n) 


1.052 
1.341 


1.118 
1.515 


1.166 
1.642 


1.206 
1.741 
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1.819 


1.264 
1.883 
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FIG. 22: Lower bounds on BH production cross sections for n = 2, ... ,7 from below, assuming 
Md = I TeV and 

Xmin — 1- Energy loss has been included according to Eq. (85). The SM 
cross section cr^^, as given in Eq. (75), is indicated by the dotted line. The typical range of 
pp interactions, as well as cross sections required for shower triggering in different characteristic 
targets are also shown for comparison. This figure is courtesy of Jonathan Feng. 



where i labels parton species and the fi{x, Q) are pdf 's (to derived the BH production cross 
sections shown in Fig. 22 we used the CTEQ6M pdf's [281, 282]). The momentum scale Q is 
taken as r~^, which is a typical momentum transfer during the gravitational collapse [283]. 
In contrast to SM processes, BH production is not supressed by perturbative couplings and 
is enhanced by the sum over all partons, particularly the gluon. 

Subsequent to formation, the BH proceeds to decay dominantly through radiation of 
standard SM particles on the 3-brane [284]. This occurs because the SM particles live on 
the brane, so that the relevant phase space for BH decay into these fields is governed by 
the 4-dimensional projection of the horizon area of the (4 + n)-dimensional BH, and by 
the Hawking temperature which is common to bulk and brane. Modulo some grey body 
factors [285-287], the dominance of the SM radiation is a result of the much larger number 
of degrees of freedom. 

The choice of a lower limit of integration in Eq. (85) requires additional explanation. 
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This limit determines the minimum mass for BHs that will be included in the calculation. 
The semiclassical description outlined above is only reliable when the energy of the emitted 
particle is small compared to the BH mass, because it is only under this condition that both 
the gravitational field of the brane and the back reaction of the metric during the emission 
process can be safely neglected [288]. Since the total number of particles emitted by the BH 
is roughly equal to its entropy, 

47rMBHr,(MBH) 

^- ^T2 ' ^^^^ 

the criterion we employ is to assume that the BH has an entropy -S" ^ 1. For Xmin = 
-^BH,min/A^r) > 3 and n > 5, one finds S > 10, so that most of the decay process can be 
well described within the semiclassical approximation. Moreover, the string cross section 
derived [289] from the Virasoro-Shapiro amplitude is expected to be considerably larger than 
that given in Eq. (82), and so it may be reasonable even to take Xmin as low as 1, for which 
^>3[290]. 

Although the BH production cross section, 0{M^^), is about 5 orders of magnitude 
smaller than QCD cross sections, (9(Aqqj-,), one of the most startling prdictions of TeV- 
scale gravity theories is that at the LHC, BH events could be filtered out of the QCD 
background, both in pp [272, 273] (see also [291-296]) as well as in PbPb [297] collisions. 
At energies of interest, however, the cosmic ray luminosity, L ~ 10^^^ cm~^ s~"^, is about 50 
orders of magnitude smaller than the LHC luminosity, thus making it futile to hunt for BHs 
in baryonic cosmic rays. On the other hand, as can be seen in Fig. 22, although greatly re- 
duced by the cross section for BH production, neutrino interaction lengths are still far larger 
than the Earth's atmospheric depth. Neutrinos therefore would produce BHs with roughly 
equal probability at any point in the atmosphere. As a result, the light descendants of the 
BH may initiate low-altitude, quasi-horizontal showers at rates significantly higher than SM 
predictions. Because of these considerations the atmosphere provides a buffer against con- 
tamination by mismeasured baryons (for which the electromagnetic channel is filtered out) 
allowing a good characterization of BH-induccd showers when S ^ 1 [298-302].^^ Further- 
more, a similar technique to that employed in discriminating between photon and hadron 
showers can be applied to isolate BH mediated showers from neutrino SM events [306]. 
Specifically, if an anomalously large quasi-horizontal deep shower rate is found, it may be 
ascribed to either an enhancement of the incoming neutrino fiux, or an enhancement in the 
neutrino-nucleon cross section. However, these two possibilities may be distinguished by 
separately binning events which arrive at very small angles to the horizontal, the so-called 
"Earth-skimming" events. An enhanced flux will increase both quasi-horizontal and Earth- 
skimming event rates, whereas a large BH cross section suppresses the latter, because the 
hadronic decay products of BH evaporation do not escape the Earth's crust. 

In summary, the signal for ultra high energy neutrinos is quasi-horizontal giant air showers 
initiated deep in the atmosphere: showers with large electromagnetic components, curved 
fronts, and signals well spread over time. These shower characteristics are easily differ- 
entiated from EAS initiated by baryons or photons. The low target density for neutrino 
interactions provided by the atmosphere must be compensated by monitoring large areas 
at the Earth's surface. In particular, the PAO will have an acceptance exceeding 1 km^ 



Additionally, neutrinos that traverse the atmosphere unscathed may produce BHs through interactions in 
the ice or water and be detected by neutrino telescopes [303-305]. 
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of water for > 10^ GeV, and thus will be able to search for extraterrestrial sources of 
ultra high energy neutrinos. Moreover, this observatory holds great promise for probing 
physics beyond the SM. An optimist might even imagine the discovery of microscopic BHs, 
the telltale signature of the universe's unseen dimensions. 

VI. EAS IN A NUTSHELL 

In this article, we have reviewed the general properties and techniques for modelling air 
showers initiated by ultra high energy particles interacting in the Earth's atmosphere. 

The incidence of a single high energy particle on the upper atmosphere gives rise to a 
roughly conical cascade of particles which reaches the Earth in the form of a giant "saucer" 
traveling at nearly the speed of light. The number of secondaries in the cascade readily 
increases through subsequent generations of particle interactions. Because of the prompt 
decay of neutral pions, about 30% of the energy in each generation is transferred to an 
electromagnetic cascade. Roughly speaking, at 10^^ GeV, baryons and charged pions have 
interaction lengths of the order of 40 g/cm^, increasing to about 60 g/cm^ at 10^ GeV. Ad- 
ditionally, below 10^° GeV, photons, electrons, and positrons have mean interaction lengths 
of 37 g/cm^, whereas above this critical energy the competing LPM and geomagnetic effects 
lead to interaction lengths between 45 — 60 g/cm^. Altogether, the atmosphere acts as a nat- 
ural colorimeter with variable density, providing a vertical thickness of 26 radiation lengths 
and about 15 interaction lengths. Amusingly, this is not too different from the number of 
radiation and interaction lengths at the LHC detectors. 

The number of muons does not increase linearly with energy, because of the previously 
mentioned pionization process: at higher energy more generations are required to cool the 
pions to the point where they are likely to decay before interaction. Production of extra 
generations results in a larger fraction of the energy being lost to the electromagnetic cascade, 
and hence a smaller fraction of the original energy being delivered to the tt^. Ultimately, 
about 90% of the primary particle's energy is dissipated in the electromagnetic cascade. 
The remaining energy is carried by hadrons, as well as muons and neutrinos produced in tt^ 
decays. 

By the time they reach the ground, relatively vertical showers have evolved fronts with a 
curvature radius of a few km, and far from the shower core their constituent particles are well 
spread over time, typically of the order of a few microseconds. For such a shower both the 
muon component and a large portion of the electromagnetic component survive to reach the 
ground, and their lateral distributions can be accurately parametrized. Although the lateral 
distribution function depends on the experiment, surface measurements of both gamma- 
and baryon-induced showers can be fitted with NKG-hke formulae. Prom this distribution 
the primary energy can be determined. 

For inclined showers the electromagnetic component is absorbed long before reaching the 
ground, as it has passed through the equivalent of several vertical atmospheres: 2 at a zenith 
angle of 60°, 3 at 70°, and 6 at 80°. In these showers, only high energy muons created in the 
first few generations of particles survive past 2 equivalent vertical atmospheres. The rate 
of energy attenuation for muons is much smaller than it is for electrons, thus the shape of 



For example, the CMS electromagnetic calorimeter is > 25 radiation lengths deep, and the hadron 
calorimeter constitutes 11 interaction lengths. 
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the resulting shower front is very flat (with curvature radius above 100 km), and its time 
extension is very short (less than 50 ns). The damping of the electromagnetic component of 
the shower provides a means to search for both photon and neutrino primaries. In particular 
quasi-horizontal showers with electromagnetic components at the ground would suggest a 
deeply penetrating primary, such as the elusive ultra high energy neutrino. 

As we have seen, the chief uncertainty in shower modelling arises from lack of definitive 
knowledge about the nature of hadronic interactions. This is because the center-of-mass 
energies involved in cosmic ray collisions are orders of magnitude beyond that achievable in 
present and foreseeable future experiments. Moreover, man-made accelerators are designed 
to probe QCD physics in the high tranvcrsc momentum region, and air shower physics is 
driven by interactions in the very forward direction. The analysis of extensive air showers 
then requires the extrapolation of hadronic interaction models more than 2 orders of mag- 
nitude in center-of-mass energy beyond the highest accelerator energies {^/s ~ 1.8 TeV) 
to date. In fact, the required extrapolation is much greater than this because air showers 
involve nuclei as well as single hadrons both as targets and projectiles. Efforts towards 
improving our understanding of soft and semi-hard processes are clearly required. 

The muon content of the shower tail is quite sensitive to the unknown details of hadronic 
physics at ultra high energies. This implies that attempts to extract composition informa- 
tion from measurements of muon content at ground level tend to be systcmatics dominated. 
There are, however, complimentary methods for uncovering the primary species which are 
less dependent on knowledge of the hadronic physics. One well-established method involves 
using fluorescence telescopes to determine the energy dependence of the depth of shower 
maximum, the so-called elongation rate, which is sensitive to the evolution of the primary 
composition with energy. Unfortunately, fluorescence detection has its own set of system- 
atic uncertainties associated with the knowledge of atmospheric properties as a function of 
time. Future hybrid experiments, such as the PAO, will record events with simultaneous 
observation of particles reaching the ground and the shower profile in the atmosphere, and 
thus provide a new arsenal of data for controlling the systematic errors. Furthermore, the 
giant aperture of the array will generate an extensive air shower sample of unprecedented 
size, ushering in a golden age of cosmic ray physics. 
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